import numpy as np
from numpy import *
import scipy
from scipy.optimize import minimize 
from scipy.integrate import nquad
import matplotlib.pyplot as plt
import matplotlib
import time
import sys
from numba import jit
from scipy.stats import ortho_group
import h5py
matplotlib.use('Agg')
np.set_printoptions(threshold=np.inf)
#########################################
#Use input 96 for agreement with exp work arXiv:2109.12631 (corresponds to D_0=0)
################################################################################################################################
# Number of points and iterations
nptsd=10
npts=3*nptsd**2-3*nptsd+1-2*nptsd-(nptsd-2)
eSqOverEpsilon0=4*np.pi*(1/137)*(1.97*10**(-4))
print('Using ',npts,'points in the reduced BZ')

#lattice spacing
a=1.42*10**(-10) 
N=1
# Dirac momentum of 1 layer of graphene
kD=4*np.pi/(3*a*sqrt(3))
# Twist angle
theta=1.53*np.pi/180

A=N*3*np.sqrt(3)/2*a**2/(4*np.sin(theta/2)**2)

# momentum scale of MBZ - ktheta is length of one side of hexagon of MBZ
ktheta=2*kD*np.sin(np.abs(theta)/2)
AMBZ=np.sin(theta/2)**2*(np.pi)**2*32/(3*np.sqrt(3)*a**2)

aMoire=3*np.sqrt(3)/2*a**2/(4*np.sin(theta/2)**2)
dscreen=40*10**(-9)
Prefactor=1/npts*AMBZ*1/(2*np.pi)**2

# alpha=w_{AB}/(vF*kD*sin(theta)) - dimensionaless constant to characterize magic angle
alpha=.414
# fermi velocity of dirac cone in MBZ
vF=1/(kD*theta)
hbar = 6.6*10**(-13)
vF=10**6*hbar
# Note energies are of order vF*kD*theta
# Potential parameters
epsilon=10
d=40*10**(-9)

# Useful matrices
Id2=np.block([[1,0],[0,1]])
pauliz=np.block([[1,0],[0,-1]])
paulix=np.block([[0,1],[1,0]])
pauliy=np.block([[0,-1j],[1j,0]])


#Lattice spacing and lattice spacing in reciprocal space and area of reduced BZdomainsize = (3*ktheta**2*np.sqrt(3)/2)
phaseStep=12
# 4/3 since ktheta has lattice spacing a-> 2a/sqrt(3)

# Number of HF iterations
iterations=13

print('starting index: ',int(sys.argv[1]))
# Construction of MBZ. Define primitive lattice vectors in MBZ
q1x=0
q1y=-ktheta

q2x=ktheta*np.sqrt(3)/2
q2y=ktheta*1/2

q3x=-ktheta*np.sqrt(3)/2
q3y=ktheta*1/2

# Bravais lattice vectors in MBZ
b1x=np.sqrt(3)/2*ktheta
b1y=3/2*ktheta
b2x=np.sqrt(3)/2*ktheta
b2y=-3/2*ktheta

# Number of shells used to construct MBZ - ie add another "shell" of Moire unit cells nshells*b_reciprocal from gamma point of zeroth MBZ
nshells=3
nactiveshells=1
radius=np.sqrt((nshells*b2x+q1x)**2+(nshells*b2y+q1y)**2)+.001*ktheta
shell1momentax=np.array([b1x,b2x,-b1x,-b2x,b1x+b2x,-b1x-b2x])
shell1momentay=np.array([b1y,b2y,-b1y,-b2y,b1y+b2y,-b1y-b2y])

shell2momentax=np.array([2*b1x,2*b2x,-2*b1x,-2*b2x,2*b1x+2*b2x,-2*b1x-2*b2x,2*b1x+b2x,-2*b1x-b2x,2*b2x+b1x,-2*b2x-b1x,b1x-b2x,b2x-b1x])
shell2momentay=np.array([2*b1y,2*b2y,-2*b1y,-2*b2y,2*b1y+2*b2y,-2*b1y-2*b2y,2*b1y+b2y,-2*b1y-b2y,2*b2y+b1y,-2*b2y-b1y,b1y-b2y,b2y-b1y])

# store kx and ky values of these points
qxvalstot=[]
qyvalstot=[]
# store layer index for each q point (layer 1 and 3 will share some)
qvalslayer=[]

for i in range(-10,10):
    for j in range(-10,10):
        if np.sqrt((i*b1x+j*b2x+q1x)**2+(i*b1y+j*b2y+q1y)**2)<radius:

            qxvalstot.append(np.round(i*b1x+j*b2x+q1x,3))
            qyvalstot.append(np.round(i*b1y+j*b2y+q1y,3))

            qvalslayer.append(2)
for i in range(-10,10):
    for j in range(-10,10):
        if np.sqrt((i*b1x+j*b2x+q1x+q2x)**2+(i*b1y+j*b2y+q1y+q2y)**2)<radius:

            qxvalstot.append(np.round(i*b1x+j*b2x+q1x+q2x,3))
            qyvalstot.append(np.round(i*b1y+j*b2y+q1y+q2y,3))

            qvalslayer.append(1)

for i in range(-10,10):
    for j in range(-10,10):
        if np.sqrt((i*b1x+j*b2x+q1x+q2x)**2+(i*b1y+j*b2y+q1y+q2y)**2)<radius:

            qxvalstot.append(np.round(i*b1x+j*b2x+q1x+q2x,3))
            qyvalstot.append(np.round(i*b1y+j*b2y+q1y+q2y,3))

            qvalslayer.append(3)
#total number of bands n
qxvalstotminus=[]
qyvalstotminus=[]
# store layer index for each q point (layer 1 and 3 will share some)
qvalslayerminus=[]

for i in range(np.size(qxvalstot)):
	qxvalstotminus.append(-qxvalstot[i])
	qyvalstotminus.append(-qyvalstot[i])
	qvalslayerminus.append(qvalslayer[i])


nbands=np.size(qxvalstot)*2
nbandsminus=nbands-2
nbandsplus=nbands+2
nactive=nbandsplus-nbandsminus

# Starting projector to be used to compute new projector
nactivehalf=int(nactive/2)
zeract=np.full((nactivehalf*4,nactivehalf*4),0.+0.j)
Chi = np.block([[np.identity(nactivehalf*4),zeract],[zeract,zeract]])

zerosupper2=np.block([[1,0,0,0,0,0,0,0],[0,1,0,0,0,0,0,0],[0,0,1,0,0,0,0,0],[0,0,0,1,0,0,0,0],[0,0,0,0,1,0,0,0],[0,0,0,0,0,1,0,0],[0,0,0,0,0,0,1,0],[0,0,0,0,0,0,0,1/2]])
zeroslower2=np.block([[1/2,0,0,0,0,0,0,0],[0,0,0,0,0,0,0,0],[0,0,0,0,0,0,0,0],[0,0,0,0,0,0,0,0],[0,0,0,0,0,0,0,0],[0,0,0,0,0,0,0,0],[0,0,0,0,0,0,0,0],[0,0,0,0,0,0,0,0]])

zerosupper4=np.block([[1,0,0,0,0,0,0,0],[0,1,0,0,0,0,0,0],[0,0,1,0,0,0,0,0],[0,0,0,1,0,0,0,0],[0,0,0,0,1,0,0,0],[0,0,0,0,0,1,0,0],[0,0,0,0,0,0,1/2,0],[0,0,0,0,0,0,0,1/2]])
zeroslower4=np.block([[1/2,0,0,0,0,0,0,0],[0,1/2,0,0,0,0,0,0],[0,0,0,0,0,0,0,0],[0,0,0,0,0,0,0,0],[0,0,0,0,0,0,0,0],[0,0,0,0,0,0,0,0],[0,0,0,0,0,0,0,0],[0,0,0,0,0,0,0,0]])

zerosupper6=np.block([[1,0,0,0,0,0,0,0],[0,1,0,0,0,0,0,0],[0,0,1,0,0,0,0,0],[0,0,0,1,0,0,0,0],[0,0,0,0,1,0,0,0],[0,0,0,0,0,1/2,0,0],[0,0,0,0,0,0,1/2,0],[0,0,0,0,0,0,0,1/2]])
zeroslower6=np.block([[1/2,0,0,0,0,0,0,0],[0,1/2,0,0,0,0,0,0],[0,0,1/2,0,0,0,0,0],[0,0,0,0,0,0,0,0],[0,0,0,0,0,0,0,0],[0,0,0,0,0,0,0,0],[0,0,0,0,0,0,0,0],[0,0,0,0,0,0,0,0]])

zerosupper8=np.block([[1,0,0,0,0,0,0,0],[0,1,0,0,0,0,0,0],[0,0,1,0,0,0,0,0],[0,0,0,1,0,0,0,0],[0,0,0,0,1/2,0,0,0],[0,0,0,0,0,1/2,0,0],[0,0,0,0,0,0,1/2,0],[0,0,0,0,0,0,0,1/2]])
zeroslower8=np.block([[1/2,0,0,0,0,0,0,0],[0,1/2,0,0,0,0,0,0],[0,0,1/2,0,0,0,0,0],[0,0,0,1/2,0,0,0,0],[0,0,0,0,0,0,0,0],[0,0,0,0,0,0,0,0],[0,0,0,0,0,0,0,0],[0,0,0,0,0,0,0,0]])

Chi2zeros = np.block([[zerosupper2,zeract],[zeract,zeroslower2]])
Chi4zeros = np.block([[zerosupper4,zeract],[zeract,zeroslower4]])
Chi6zeros = np.block([[zerosupper6,zeract],[zeract,zeroslower6]])
Chi8zeros = np.block([[zerosupper8,zeract],[zeract,zeroslower8]])


# Symmetries -- projected
# Number of shells used to construct MBZ - ie add another "shell" of Moire unit cells nshells*b_reciprocal from gamma point of zeroth MBZ
nshells=3
radius=np.sqrt((nshells*b2x+q1x)**2+(nshells*b2y+q1y)**2)+.001*ktheta

shell1momentax=[]
shell1momentay=[]

shell2momentax=[]
shell2momentay=[]

shell3momentax=[]
shell3momentay=[]

shell4momentax=[]
shell4momentay=[]

shell5momentax=[]
shell5momentay=[]

for a in range(-10,10):
	for b in range(-10,10):
		if np.sqrt((a*b1x+b*b2x)**2+(a*b1y+b*b2y)**2)<np.sqrt((1*b2x)**2+(1*b2y)**2)+.001*ktheta:
			if a!=0 or b!=0:

				shell1momentax.append(a*b1x+b*b2x)
				shell1momentay.append(a*b1y+b*b2y)
for a in range(-10,10):
	for b in range(-10,10):
		if np.sqrt((a*b1x+b*b2x)**2+(a*b1y+b*b2y)**2)<np.sqrt((2*b2x)**2+(2*b2y)**2)+.001*ktheta and np.sqrt((a*b1x+b*b2x)**2+(a*b1y+b*b2y)**2)>np.sqrt((1*b2x)**2+(1*b2y)**2)+.001*ktheta:
			

			shell2momentax.append(a*b1x+b*b2x)
			shell2momentay.append(a*b1y+b*b2y)

for a in range(-10,10):
	for b in range(-10,10):
		if np.sqrt((a*b1x+b*b2x)**2+(a*b1y+b*b2y)**2)<np.sqrt((3*b2x)**2+(3*b2y)**2)+.01*ktheta and np.sqrt((a*b1x+b*b2x)**2+(a*b1y+b*b2y)**2)>np.sqrt((2*b2x)**2+(2*b2y)**2)+.001*ktheta:
			
			shell3momentax.append(a*b1x+b*b2x)
			shell3momentay.append(a*b1y+b*b2y)

for a in range(-10,10):
	for b in range(-10,10):
		if np.sqrt((a*b1x+b*b2x)**2+(a*b1y+b*b2y)**2)<np.sqrt((4*b2x)**2+(4*b2y)**2)+.001*ktheta and np.sqrt((a*b1x+b*b2x)**2+(a*b1y+b*b2y)**2)>np.sqrt((3*b2x)**2+(3*b2y)**2)+.001*ktheta:
			shell4momentax.append(a*b1x+b*b2x)
			shell4momentay.append(a*b1y+b*b2y)
for a in range(-10,10):
	for b in range(-10,10):
		if np.sqrt((a*b1x+b*b2x)**2+(a*b1y+b*b2y)**2)<np.sqrt((5*b2x)**2+(5*b2y)**2)+.001*ktheta and np.sqrt((a*b1x+b*b2x)**2+(a*b1y+b*b2y)**2)>np.sqrt((4*b2x)**2+(4*b2y)**2)+.001*ktheta:
			shell5momentax.append(a*b1x+b*b2x)
			shell5momentay.append(a*b1y+b*b2y)
shell1momentaxround=np.round(shell1momentax/ktheta,3)
shell2momentaxround=np.round(shell2momentax/ktheta,3)
shell3momentaxround=np.round(shell3momentax/ktheta,3)

shell1size=np.size(shell1momentax)
shell2size=np.size(shell2momentax)
shell3size=np.size(shell3momentax)
shell1momentayround=np.round(shell1momentay/ktheta,3)
shell2momentayround=np.round(shell2momentay/ktheta,3)
shell3momentayround=np.round(shell3momentay/ktheta,3)

shell1momentax=np.array(shell1momentax)
shell2momentax=np.array(shell2momentax)
shell3momentax=np.array(shell3momentax)
shell4momentax=np.array(shell4momentax)
shell5momentax=np.array(shell5momentax)


shell1momentay=np.array(shell1momentay)
shell2momentay=np.array(shell2momentay)
shell3momentay=np.array(shell3momentay)
shell4momentay=np.array(shell4momentay)
shell5momentay=np.array(shell5momentay)

# Rotate k vectors
def C3rotate(kx,ky):
    kxprime=np.cos(2*np.pi/3)*kx-np.sin(2*np.pi/3)*ky
    kyprime=np.sin(2*np.pi/3)*kx+np.cos(2*np.pi/3)*ky
    return kxprime,kyprime
def C2zrotate(kx,ky):
    kxprime=-kx
    kyprime=-ky
    return kxprime,kyprime
def C2xrotate(kx,ky):
    kxprime=kx
    kyprime=-ky
    return kxprime,kyprime
def C2yrotate(kx,ky):
    kxprime=-kx
    kyprime=ky
    return kxprime,kyprime
# Define unitary parts of symms
def C3():
	nbands=np.size(qxvalstot)*2
	UC3=np.full((nbands,nbands),0.+0.j)
	UC3minus=np.full((nbands,nbands),0.+0.j)
	for i in range(np.size(qxvalstot)):
		kxprime,kyprime=C3rotate(qxvalstot[i],qyvalstot[i])


		x=np.where(np.round(qxvalstot/ktheta,3)==round(kxprime/ktheta, 3))
		y=np.where(np.round(qyvalstot/ktheta,3)==round(kyprime/ktheta, 3))


		for j in range( np.size(np.intersect1d(x, y))):
			if qvalslayer[i]==qvalslayer[np.intersect1d(x, y)[j]]:
				c3index=np.intersect1d(x, y)[j]

		UC3[2*c3index][2*i]=np.exp(-1j*2*np.pi/3)
		UC3[2*c3index+1][2*i+1]=np.exp(1j*2*np.pi/3)

	for i in range(np.size(qxvalstotminus)):
		kxprime,kyprime=C3rotate(qxvalstotminus[i],qyvalstotminus[i])


		x=np.where(np.round(qxvalstotminus/ktheta,3)==round(kxprime/ktheta, 3))
		y=np.where(np.round(qyvalstotminus/ktheta,3)==round(kyprime/ktheta, 3))


		for j in range( np.size(np.intersect1d(x, y))):
			if qvalslayerminus[i]==qvalslayerminus[np.intersect1d(x, y)[j]]:
				c3index=np.intersect1d(x, y)[j]

		UC3minus[2*c3index][2*i]=np.exp(1j*2*np.pi/3)
		UC3minus[2*c3index+1][2*i+1]=np.exp(-1j*2*np.pi/3)

	blockzeros=np.full((nbands,nbands),0)
	UC3Full=np.block([[UC3,blockzeros,blockzeros,blockzeros],[blockzeros,UC3minus,blockzeros,blockzeros],[blockzeros,blockzeros,UC3,blockzeros],[blockzeros,blockzeros,blockzeros,UC3minus]])



	return UC3Full
def C2z():
	nbands=np.size(qxvalstot)*2
	UC2z=np.full((nbands*2,nbands*2),0.+0.j)
	for i in range(np.size(qxvalstot)):
		kxprime,kyprime=C2zrotate(qxvalstot[i],qyvalstot[i])


		x=np.where(np.round(qxvalstotminus/ktheta,3)==round(kxprime/ktheta, 3))
		y=np.where(np.round(qyvalstotminus/ktheta,3)==round(kyprime/ktheta, 3))

		c2zindex=0

		intersect=np.intersect1d(x, y)

		for j in range( np.size(intersect)):
			if qvalslayer[i]==qvalslayerminus[np.intersect1d(x, y)[j]]:
				c2zindex=np.intersect1d(x, y)[j]


		UC2z[2*i][nbands+2*c2zindex+1]=1
		UC2z[2*i+1][nbands+2*c2zindex]=1
		UC2z[nbands+2*c2zindex][2*i+1]=1
		UC2z[nbands+2*c2zindex+1][2*i]=1



	blockzeros=np.full((2*nbands,2*nbands),0)
	UC2zFull=np.block([[UC2z,blockzeros],[blockzeros,UC2z]])
	return UC2zFull
def C2x():
    nbands=np.size(qxvalstot)*2
    UC2x=np.full((nbands,nbands),0.+0.j)
    for i in range(np.size(qxvalstot)):
        kxprime,kyprime=C2xrotate(qxvalstot[i],qyvalstot[i])


        x=np.where(np.round(qxvalstot/ktheta,3)==round(kxprime/ktheta, 3))
        y=np.where(np.round(qyvalstot/ktheta,3)==round(kyprime/ktheta, 3))

        c2zindex=0

        for j in range( np.size(np.intersect1d(x, y))):
            if qvalslayer[i]==qvalslayer[np.intersect1d(x, y)[j]]:
                c2zindex=np.intersect1d(x, y)[j]

        UC2x[2*i][2*i+1]=1
        UC2x[2*i+1][2*i]=1


    blockzeros=np.full((nbands,nbands),0)
    UC2xFull=np.block([[UC2x,blockzeros,blockzeros,blockzeros],[blockzeros,UC2x,blockzeros,blockzeros],[blockzeros,blockzeros,UC2x,blockzeros],[blockzeros,blockzeros,blockzeros,UC2x]])
    return UC2xFull
def Chiral():
	base=np.full((npts,nactive*4,nactive*4),0.+0.j)

	blockivc=np.full((nbands,nbands),0.+0.j)
	for i in range(np.size(qxvalstot)):
		blockivc[2*i][2*i]=1
		blockivc[2*i+1][2*i+1]=-1
	zerblock=np.full((nbands,nbands),0.+0.j)
	UnProjected=np.block([[blockivc,zerblock,zerblock,zerblock],[zerblock,blockivc,zerblock,zerblock],[zerblock,zerblock,blockivc,zerblock],[zerblock,zerblock,zerblock,blockivc]])

	return UnProjected
def T():
    nbands=np.size(qxvalstot)*2
    UC2z=np.full((nbands,nbands),0.+0.j)
    for i in range(np.size(qxvalstot)):
        kxprime,kyprime=C2zrotate(qxvalstot[i],qyvalstot[i])


        x=np.where(np.round(qxvalstot/ktheta,3)==round(kxprime/ktheta, 3))
        y=np.where(np.round(qyvalstot/ktheta,3)==round(kyprime/ktheta, 3))

        c2zindex=0

        for j in range( np.size(np.intersect1d(x, y))):
            if qvalslayer[i]==qvalslayer[np.intersect1d(x, y)[j]]:
                c2zindex=np.intersect1d(x, y)[j]

        UC2z[2*i][2*i]=1
        UC2z[2*i+1][2*i+1]=1


    blockzeros=np.full((nbands,nbands),0)
    UC2zFull=np.block([[blockzeros,UC2z,blockzeros,blockzeros],[UC2z,blockzeros,blockzeros,blockzeros],[blockzeros,blockzeros,blockzeros,UC2z],[blockzeros,blockzeros,UC2z,blockzeros]])
    return UC2zFull
def T0():
	nbands=np.size(qxvalstot)*2
	UT=np.full((nbands*2,nbands*2),0.+0.j)
	for i in range(np.size(qxvalstot)):
		kxprime,kyprime=C2zrotate(qxvalstot[i],qyvalstot[i])


		x=np.where(np.round(qxvalstotminus/ktheta,3)==round(kxprime/ktheta, 3))
		y=np.where(np.round(qyvalstotminus/ktheta,3)==round(kyprime/ktheta, 3))

		c2zindex=0

		intersect=np.intersect1d(x, y)

		for j in range( np.size(intersect)):
			if qvalslayer[i]==qvalslayerminus[np.intersect1d(x, y)[j]]:
				c2zindex=np.intersect1d(x, y)[j]


		UT[2*i][nbands+2*c2zindex]=1
		UT[2*i+1][nbands+2*c2zindex+1]=1
		UT[nbands+2*c2zindex][2*i]=1
		UT[nbands+2*c2zindex+1][2*i+1]=1


	blockzeros=np.full((2*nbands,2*nbands),0)
	return np.block([[UT,blockzeros],[blockzeros,UT]])
def T0Spin():
	nbands=np.size(qxvalstot)*2
	UT=np.full((nbands*2,nbands*2),0.+0.j)
	for i in range(np.size(qxvalstot)):
		kxprime,kyprime=C2zrotate(qxvalstot[i],qyvalstot[i])


		x=np.where(np.round(qxvalstotminus/ktheta,3)==round(kxprime/ktheta, 3))
		y=np.where(np.round(qyvalstotminus/ktheta,3)==round(kyprime/ktheta, 3))

		c2zindex=0

		intersect=np.intersect1d(x, y)

		for j in range( np.size(intersect)):
			if qvalslayer[i]==qvalslayerminus[np.intersect1d(x, y)[j]]:
				c2zindex=np.intersect1d(x, y)[j]


		UT[2*i][nbands+2*c2zindex]=1
		UT[2*i+1][nbands+2*c2zindex+1]=1
		UT[nbands+2*c2zindex][2*i]=1
		UT[nbands+2*c2zindex+1][2*i+1]=1


	blockzeros=np.full((2*nbands,2*nbands),0)
	return 1j*np.block([[blockzeros,-1j*UT],[1j*UT,blockzeros]])
def U1v():
    nbands=np.size(qxvalstot)*2
    UU1v=np.full((nbands,nbands),0.+0.j)
    for i in range(np.size(qxvalstot)):
        UU1v[2*i][2*i]=1
        UU1v[2*i+1][2*i+1]=1


    blockzeros=np.full((nbands,nbands),0)
    UU1vFull=np.block([[UU1v,blockzeros,blockzeros,blockzeros],[blockzeros,-UU1v,blockzeros,blockzeros],[blockzeros,blockzeros,UU1v,blockzeros],[blockzeros,blockzeros,blockzeros,-UU1v]])
    return UU1vFull
def Spin():
    nbands=np.size(qxvalstot)*2
    Uspin=np.full((nbands,nbands),0.+0.j)
    for i in range(np.size(qxvalstot)):
        Uspin[2*i][2*i]=1
        Uspin[2*i+1][2*i+1]=1


    blockzeros=np.full((nbands,nbands),0)
    UspinFull=np.block([[Uspin,blockzeros,blockzeros,blockzeros],[blockzeros,Uspin,blockzeros,blockzeros],[blockzeros,blockzeros,-Uspin,blockzeros],[blockzeros,blockzeros,blockzeros,-Uspin]])
    return UspinFull
def ChiralSymmetry(P):
	total=0

	for i in range(npts):

		Pcheck=P[i]


		total=total+np.sum(np.abs(np.round(Chiralmatrix.dot(P[i]).dot(Chiralmatrix)-Pcheck,6)))


		print('testing: ',i,np.sum(np.abs(np.round((Chiralmatrix.dot(P[i]).dot(Chiralmatrix)-Pcheck),6))))


	if total<10**(-3):
	    return 1
	else:
	    return 0
# Symmetry action on projector
def C3Symmetry(P):
	total=0

	for i in range(npts):

		Pcheck=P[kc3index[i]]


		# total=total+np.sum(np.abs(np.round(np.conjugate(C3matrix).dot(P[i]).dot(np.transpose(C3matrix))-Pcheck,6)))


		print('testing: ',i,np.sum(np.abs(np.round((np.conjugate(C3matrix).dot(P[i]).dot(np.transpose(C3matrix))-Pcheck),20))))


	if total<10**(-3):
	    return 1
	else:
	    return 0
def C2zSymmetry(P):
	total=0

	for i in range(npts):
		if np.abs(kxvals[i]+b1x)<10**(-3)*ktheta or np.abs(kyvals[i]-kxvals[i]*q2y/q2x-ktheta)<10**(-3) or np.abs(kyvals[i]+kxvals[i]*q2y/q2x-ktheta)<10**(-3):
			continue
		Pcheck=P[kminusindex[i]]


		# total=total+np.sum(np.abs(np.round(np.conjugate(C2zmatrix).dot(P[i]).dot(np.transpose(C2zmatrix))-Pcheck,6)))


		print('testing: ',i,np.sum(np.abs(np.round(np.conjugate(C2zmatrix).dot(P[i]).dot(np.transpose(C2zmatrix))-Pcheck,6))))


	if total<10**(-3):
	    return 1
	else:
	    return 0
def C2zTSymmetry(P):
	total=0

	for i in range(npts):
		if np.abs(kxvals[i]+b1x)<10**(-3)*ktheta or np.abs(kyvals[i]-kxvals[i]*q2y/q2x-ktheta)<10**(-3) or np.abs(kyvals[i]+kxvals[i]*q2y/q2x-ktheta)<10**(-3):
			continue
		Pcheck=P[kminusindex[i]]


		total=total+np.sum(np.abs(np.round(np.conjugate(C2zmatrix.dot(T0matrix)).dot(P[i]).dot(np.transpose(C2zmatrix.dot(T0matrix)))-Pcheck,6)))


		print('testing: ',i,np.sum(np.abs(np.round(np.conjugate(C2zmatrix.dot(T0matrix)).dot(np.conjugate(P[i])).dot(np.transpose(C2zmatrix.dot(T0matrix)))-Pcheck,6))))


	if total<10**(-3):
	    return 1
	else:
	    return 0
def C2zTSymmetry(P):
	total=0

	for i in range(npts):

		Pcheck=P[i]


		total=total+np.sum(np.abs(np.round(np.conjugate(C2zmatrix).dot(P[i]).dot(np.transpose(C2zmatrix))-Pcheck,6)))


		print('testing: ',i,np.sum(np.abs(np.round(np.conjugate(C2zmatrix.dot(T0matrix)).dot(np.conjugate(P[i])).dot(np.transpose(C2zmatrix.dot(T0matrix)))-Pcheck,6))))


	if total<10**(-3):
	    return 1
	else:
	    return 0
def TSymmetry(P):
	total=0

	for i in range(npts):
		if np.abs(kxvals[i]+b1x)<10**(-3)*ktheta or np.abs(kyvals[i]-kxvals[i]*q2y/q2x-ktheta)<10**(-3) or np.abs(kyvals[i]+kxvals[i]*q2y/q2x-ktheta)<10**(-3):
			continue
		Pcheck=P[kminusindex[i]]




		print('testing: ',i,np.sum(np.abs(np.round(np.conjugate(T0matrix).dot(np.conjugate(P[i])).dot(np.transpose(T0matrix))-Pcheck,6))))


	if total<10**(-3):
	    return 1
	else:
	    return 0
def TSpinSymmetry(P):
	total=0

	for i in range(npts):
		if np.abs(kxvals[i]+b1x)<10**(-3)*ktheta or np.abs(kyvals[i]-kxvals[i]*q2y/q2x-ktheta)<10**(-3) or np.abs(kyvals[i]+kxvals[i]*q2y/q2x-ktheta)<10**(-3):
			continue
		Pcheck=P[kminusindex[i]]




		print('testing: ',i,np.sum(np.abs(np.round(np.conjugate(T0Spinmatrix).dot(np.conjugate(P[i])).dot(np.transpose(T0Spinmatrix))-Pcheck,16))))


	if total<10**(-3):
	    return 1
	else:
	    return 0
def U1vSymmetry(P):
	total=0

	for i in range(npts):
		if np.abs(kxvals[i]+b1x)<10**(-3)*ktheta or np.abs(kyvals[i]-kxvals[i]*q2y/q2x-ktheta)<10**(-3) or np.abs(kyvals[i]+kxvals[i]*q2y/q2x-ktheta)<10**(-3):
			continue
		Pcheck=P[i]


		total=total+np.sum(np.abs(np.round(np.conjugate(U1vmatrix).dot(P[i]).dot(np.transpose(U1vmatrix))-Pcheck,6)))


		print('testing: ',i,np.sum(np.abs(np.round(np.conjugate(U1vmatrix).dot(P[i]).dot(np.transpose(U1vmatrix))-Pcheck,6))))


	if total<10**(-3):
	    return 1
	else:
	    return 0
def SpinSymmetry(P):
	total=0

	for i in range(npts):
		if np.abs(kxvals[i]+b1x)<10**(-3)*ktheta or np.abs(kyvals[i]-kxvals[i]*q2y/q2x-ktheta)<10**(-3) or np.abs(kyvals[i]+kxvals[i]*q2y/q2x-ktheta)<10**(-3):
			continue
		Pcheck=P[i]


		total=total+np.sum(np.abs(np.round(np.conjugate(Spinmatrix).dot(P[i]).dot(np.transpose(Spinmatrix))-Pcheck,6)))


		print('testing: ',i,np.sum(np.abs(np.round(np.conjugate(Spinmatrix).dot(P[i]).dot(np.transpose(Spinmatrix))-Pcheck,6))))


	if total<10**(-3):
	    return 1
	else:
	    return 0
# Symmetry action on Hamiltonian
def C2zSymmetryH(P):
	total=0

	for i in range(npts):
		if np.abs(kxvals[i]+b1x)<10**(-3)*ktheta or np.abs(kyvals[i]-kxvals[i]*q2y/q2x-ktheta)<10**(-3) or np.abs(kyvals[i]+kxvals[i]*q2y/q2x-ktheta)<10**(-3) or i%3!=0:
			continue
		Pcheck=P[kminusindex[i]]




		print('testing: ',i,np.sum(np.abs(np.round(C2zmatrix.dot(P[i]).dot(np.transpose(np.conjugate(C2zmatrix)))-Pcheck,12))))


	if total<10**(-3):
	    return 1
	else:
	    return 0
def TSymmetryH(P):
	total=0

	for i in range(npts):
		if i%3!=0:
			continue

		Pcheck=P[kminusindex[i]]




		print('testing: ',i,np.sum(np.abs(np.round((T0matrix.dot(np.conjugate(P[i])).dot(np.conjugate(np.transpose(T0matrix)))-Pcheck),12))))


	if total<10**(-3):
	    return 1
	else:
	    return 0
def TSpinSymmetryH(P):
	total=0

	for i in range(npts):

		Pcheck=P[kminusindex[i]]




		print('testing: ',i,np.sum(np.abs(np.round((T0Spinmatrix.dot(np.conjugate(P[i])).dot(np.conjugate(np.transpose(T0Spinmatrix)))-Pcheck),12))))


	if total<10**(-3):
	    return 1
	else:
	    return 0
def C3SymmetryH(P):
	total=0

	for i in range(npts):
		if i%3!=0:
			continue

		Pcheck=P[kc3index[i]]




		print('testing: ',i,np.sum(np.abs(np.round((C3matrix.dot(P[i]).dot(np.conjugate(np.transpose(C3matrix)))-Pcheck),12))))


	if total<10**(-3):
	    return 1
	else:
	    return 0
# Orders
def C3Spin():
	nbands=np.size(qxvalstot)*2
	UC3=np.full((nbands,nbands),0.+0.j)
	UC3minus=np.full((nbands,nbands),0.+0.j)
	for i in range(np.size(qxvalstot)):
		kxprime,kyprime=C3rotate(qxvalstot[i],qyvalstot[i])


		x=np.where(np.round(qxvalstot/ktheta,3)==round(kxprime/ktheta, 3))
		y=np.where(np.round(qyvalstot/ktheta,3)==round(kyprime/ktheta, 3))


		for j in range( np.size(np.intersect1d(x, y))):
			if qvalslayer[i]==qvalslayer[np.intersect1d(x, y)[j]]:
				c3index=np.intersect1d(x, y)[j]

		UC3[2*c3index][2*i]=np.exp(-1j*2*np.pi/3)
		UC3[2*c3index+1][2*i+1]=np.exp(1j*2*np.pi/3)

	for i in range(np.size(qxvalstotminus)):
		kxprime,kyprime=C3rotate(qxvalstotminus[i],qyvalstotminus[i])


		x=np.where(np.round(qxvalstotminus/ktheta,3)==round(kxprime/ktheta, 3))
		y=np.where(np.round(qyvalstotminus/ktheta,3)==round(kyprime/ktheta, 3))


		for j in range( np.size(np.intersect1d(x, y))):
			if qvalslayerminus[i]==qvalslayerminus[np.intersect1d(x, y)[j]]:
				c3index=np.intersect1d(x, y)[j]

		UC3minus[2*c3index][2*i]=np.exp(1j*2*np.pi/3)
		UC3minus[2*c3index+1][2*i+1]=np.exp(-1j*2*np.pi/3)

	blockzeros=np.full((nbands,nbands),0)
	UC3Full=np.block([[UC3,blockzeros,blockzeros,blockzeros],[blockzeros,UC3minus,blockzeros,blockzeros],[blockzeros,blockzeros,UC3,blockzeros],[blockzeros,blockzeros,blockzeros,UC3minus]])



	return UC3Full
def C2zSpin():
	nbands=np.size(qxvalstot)*2
	UC2z=np.full((nbands*2,nbands*2),0.+0.j)
	for i in range(np.size(qxvalstot)):
		kxprime,kyprime=C2zrotate(qxvalstot[i],qyvalstot[i])


		x=np.where(np.round(qxvalstotminus/ktheta,3)==round(kxprime/ktheta, 3))
		y=np.where(np.round(qyvalstotminus/ktheta,3)==round(kyprime/ktheta, 3))

		c2zindex=0

		intersect=np.intersect1d(x, y)

		for j in range( np.size(intersect)):
			if qvalslayer[i]==qvalslayerminus[np.intersect1d(x, y)[j]]:
				c2zindex=np.intersect1d(x, y)[j]


		UC2z[2*i][nbands+2*c2zindex+1]=1
		UC2z[2*i+1][nbands+2*c2zindex]=1
		UC2z[nbands+2*c2zindex][2*i+1]=1
		UC2z[nbands+2*c2zindex+1][2*i]=1





	blockzeros=np.full((2*nbands,2*nbands),0)
	UC2zFull=np.block([[UC2z,blockzeros],[blockzeros,UC2z]])
	return UC2zFull


def Identity():
	base=np.full((npts,nactive*4,nactive*4),0.+0.j)

	blockspin=np.identity(nbands*2)
	zerblock=np.full((nbands*2,nbands*2),0.+0.j)
	UnProjected=np.block([[blockspin,zerblock],[zerblock,blockspin]])
	Identity=np.identity(4*nbands)
	IdblockSpecial=np.block([[1,0,0,0],[0,1,0,0],[0,0,1,0],[0,0,0,1]])
	FilledBottomBands=np.block([[2,0,0,0],[0,2,0,0],[0,0,0,0],[0,0,0,0]])
	Idblock=np.block([[2,0,0,0],[0,1,0,0],[0,0,1,0],[0,0,0,0]])
	zersmall=np.full((nactive,nactive),0.+0.j)
	for i in range(npts):

		U=HtotalEigVecs[i]

		base[i]=-1/2*np.identity(nactive*4)


	return base
def LowerBands():
	base=np.full((npts,nactive*4,nactive*4),0.+0.j)

	blockspin=np.identity(nbands*2)
	zerblock=np.full((nbands*2,nbands*2),0.+0.j)
	UnProjected=np.block([[blockspin,zerblock],[zerblock,blockspin]])
	Identity=np.identity(4*nbands)
	IdblockSpecial=np.block([[1,0,0,0],[0,1,0,0],[0,0,1,0],[0,0,0,1]])
	FilledBottomBands=np.block([[2,0,0,0],[0,2,0,0],[0,0,0,0],[0,0,0,0]])
	Idblock=np.block([[2,0,0,0],[0,1,0,0],[0,0,1,0],[0,0,0,0]])
	zersmall=np.full((nactive,nactive),0.+0.j)
	for i in range(npts):

		U=HtotalEigVecs[i]

		base[i]=np.identity(nactive*4)
		for l in range(4*nactive):
		    if l%4==2 or l%4==3:
		        base[i][l][l]=0

		
		if i==0 or i==nptsd-1:
			base[i]=1/2*np.identity(nactive*4)

	return base
def Identityshell1():
	base=np.full((np.size(shell1momentax),npts,nactive*4,nactive*4),0.+0.j)

	blockspin=np.identity(nbands*2)
	zerblock=np.full((nbands*2,nbands*2),0.+0.j)
	UnProjected=np.block([[blockspin,zerblock],[zerblock,blockspin]])
	Identity=np.identity(4*nbands)
	IdblockSpecial=np.block([[1,0,0,0],[0,1,0,0],[0,0,1,0],[0,0,0,1]])
	Idblock=np.block([[2,0,0,0],[0,1,0,0],[0,0,1,0],[0,0,0,0]])
	FilledBottomBands=np.block([[2,0,0,0],[0,2,0,0],[0,0,0,0],[0,0,0,0]])
	zersmall=np.full((nactive,nactive),0.+0.j)
	for i in range(npts):
		for k in range(np.size(shell1momentax)):

			U=HtotalEigVecsShell1[k][i]
			base[k][i]=-1/2*np.identity(nactive*4)


	return base
def Identityshell2():
	base=np.full((np.size(shell2momentax),npts,nactive*4,nactive*4),0.+0.j)

	blockspin=np.identity(nbands*2)
	zerblock=np.full((nbands*2,nbands*2),0.+0.j)
	UnProjected=np.block([[blockspin,zerblock],[zerblock,blockspin]])
	Identity=np.identity(4*nbands)
	IdblockSpecial=np.block([[1,0,0,0],[0,1,0,0],[0,0,1,0],[0,0,0,1]])
	Idblock=np.block([[2,0,0,0],[0,1,0,0],[0,0,1,0],[0,0,0,0]])
	FilledBottomBands=np.block([[2,0,0,0],[0,2,0,0],[0,0,0,0],[0,0,0,0]])
	zersmall=np.full((nactive,nactive),0.+0.j)
	for i in range(npts):
		for k in range(np.size(shell2momentax)):

		    U=HtotalEigVecsShell2[k][i]
		    base[k][i]=-1/2*np.identity(nactive*4)



	return base
# Fix order edges
def MakeShell1P(Proj):
	base=np.full((np.size(shell1momentax),npts,nactive*4,nactive*4),0.+0.j)

	for i in range(np.size(shell1momentax)):

		for j in range(npts):
			if not(np.abs(kxvals[j]+b1x)<10**(-3)*ktheta or np.abs(kyvals[j]-kxvals[j]*q2y/q2x-ktheta)<10**(-3) or np.abs(kyvals[j]+kxvals[j]*q2y/q2x-ktheta)<10**(-3)):
				continue
			U=HtotalEigVecsShell1[i][j]


			

			if (kxvals[j]+shell1momentax[i])/ktheta<0-10**(-3) and (kxvals[j]+shell1momentax[i])/ktheta*q3y/q3x-(kyvals[j]+shell1momentay[i])/ktheta<10**(-3):



				base[i][j]=np.transpose(U).dot(np.conjugate(HtotalEigVecsShell1[i][j]).dot(Proj[j]).dot(np.transpose(HtotalEigVecsShell1[i][j]))).dot(np.conjugate(U))
				

				


			elif (kxvals[j]+shell1momentax[i])/ktheta>0-10**(-3) and -(kxvals[j]+shell1momentax[i])/ktheta*q3y/q3x-(kyvals[j]+shell1momentay[i])/ktheta<-10**(-3):



				kxprime,kyprime=C3rotate((kxvals[j]+shell1momentax[i]),(kyvals[j]+shell1momentay[i]))
				kxprime,kyprime=C3rotate(kxprime,kyprime)
				kxprime,kyprime=C2zrotate(kxprime,kyprime)

				kxfind=np.where(np.round(kxvals/ktheta,3)==np.round((kxprime)/ktheta,3))
				kyfind=np.where(np.round(kyvals/ktheta,3)==np.round((kyprime)/ktheta,3))
				indexfind=np.intersect1d(kxfind,kyfind)

				if np.size(indexfind)==0:
					for k in range(np.size(shell1momentax)):
						kxfind=np.where(np.round(kxvals/ktheta,3)==np.round((kxprime-shell1momentax[k])/ktheta,3))
						kyfind=np.where(np.round(kyvals/ktheta,3)==np.round((kyprime-shell1momentay[k])/ktheta,3))
						indexfind=np.intersect1d(kxfind,kyfind)
						if np.size(indexfind)>0:
							indexnew=indexfind[0]
							kxprime=kxvals[indexnew]+shell1momentax[k]
							kyprime=kyvals[indexnew]+shell1momentay[k]
							Pmatch=np.conjugate(HtotalEigVecsShell1[k][indexnew]).dot(Proj[indexnew]).dot(np.transpose(HtotalEigVecsShell1[k][indexnew]))
					for k in range(np.size(shell2momentax)):
						kxfind=np.where(np.round(kxvals/ktheta,3)==np.round((kxprime-shell2momentax[k])/ktheta,3))
						kyfind=np.where(np.round(kyvals/ktheta,3)==np.round((kyprime-shell2momentay[k])/ktheta,3))
						indexfind=np.intersect1d(kxfind,kyfind)
						if np.size(indexfind)>0:
							indexnew=indexfind[0]
							kxprime=kxvals[indexnew]+shell2momentax[k]
							kyprime=kyvals[indexnew]+shell2momentay[k]
							Pmatch=np.conjugate(HtotalEigVecsShell2[k][indexnew]).dot(Proj[indexnew]).dot(np.transpose(HtotalEigVecsShell2[k][indexnew]))
				else:
					indexnew=indexfind[0]
					kxprime=kxvals[indexnew]
					kyprime=kyvals[indexnew]
					Pmatch=np.conjugate(HtotalEigVecs[indexnew]).dot(Proj[indexnew]).dot(np.transpose(HtotalEigVecs[indexnew]))



				base[i][j]=np.transpose(U).dot(np.conjugate(C3C2product).dot(Pmatch).dot(np.transpose(C3C2product))).dot(np.conjugate(U))




			elif (kxvals[j]+shell1momentax[i])/ktheta>0+10**(-3) and -(kxvals[j]+shell1momentax[i])/ktheta*q3y/q3x-(kyvals[j]+shell1momentay[i])/ktheta>-10**(-3) and (kxvals[j]+shell1momentax[i])/ktheta*q3y/q3x-(kyvals[j]+shell1momentay[i])/ktheta<-10**(-3):

				kxprime,kyprime=C3rotate((kxvals[j]+shell1momentax[i]),(kyvals[j]+shell1momentay[i]))

				kxfind=np.where(np.round(kxvals/ktheta,3)==np.round((kxprime)/ktheta,3))
				kyfind=np.where(np.round(kyvals/ktheta,3)==np.round((kyprime)/ktheta,3))
				indexfind=np.intersect1d(kxfind,kyfind)

				if np.size(indexfind)==0:
					for k in range(np.size(shell1momentax)):
						kxfind=np.where(np.round(kxvals/ktheta,3)==np.round((kxprime-shell1momentax[k])/ktheta,3))
						kyfind=np.where(np.round(kyvals/ktheta,3)==np.round((kyprime-shell1momentay[k])/ktheta,3))
						indexfind=np.intersect1d(kxfind,kyfind)
						if np.size(indexfind)>0:
							indexnew=indexfind[0]
							kxprime=kxvals[indexnew]+shell1momentax[k]
							kyprime=kyvals[indexnew]+shell1momentay[k]
							Pmatch=np.conjugate(HtotalEigVecsShell1[k][indexnew]).dot(Proj[indexnew]).dot(np.transpose(HtotalEigVecsShell1[k][indexnew]))
					for k in range(np.size(shell2momentax)):
						kxfind=np.where(np.round(kxvals/ktheta,3)==np.round((kxprime-shell2momentax[k])/ktheta,3))
						kyfind=np.where(np.round(kyvals/ktheta,3)==np.round((kyprime-shell2momentay[k])/ktheta,3))
						indexfind=np.intersect1d(kxfind,kyfind)
						if np.size(indexfind)>0:
							indexnew=indexfind[0]
							kxprime=kxvals[indexnew]+shell2momentax[k]
							kyprime=kyvals[indexnew]+shell2momentay[k]
							Pmatch=np.conjugate(HtotalEigVecsShell2[k][indexnew]).dot(Proj[indexnew]).dot(np.transpose(HtotalEigVecsShell2[k][indexnew]))
				else:
					indexnew=indexfind[0]
					kxprime=kxvals[indexnew]
					kyprime=kyvals[indexnew]
					Pmatch=np.conjugate(HtotalEigVecs[indexnew]).dot(Proj[indexnew]).dot(np.transpose(HtotalEigVecs[indexnew]))

				base[i][j]=np.transpose(U).dot(np.conjugate(C3C3product).dot(Pmatch).dot(np.transpose(C3C3product))).dot(np.conjugate(U))

			elif (kxvals[j]+shell1momentax[i])/ktheta>0+10**(-3) and (kxvals[j]+shell1momentax[i])/ktheta*q3y/q3x-(kyvals[j]+shell1momentay[i])/ktheta>-10**(-3):

				kxprime,kyprime=C2zrotate((kxvals[j]+shell1momentax[i]),(kyvals[j]+shell1momentay[i]))

				kxfind=np.where(np.round(kxvals/ktheta,3)==np.round((kxprime)/ktheta,3))
				kyfind=np.where(np.round(kyvals/ktheta,3)==np.round((kyprime)/ktheta,3))
				indexfind=np.intersect1d(kxfind,kyfind)

				if np.size(indexfind)==0:
					for k in range(np.size(shell1momentax)):
						kxfind=np.where(np.round(kxvals/ktheta,3)==np.round((kxprime-shell1momentax[k])/ktheta,3))
						kyfind=np.where(np.round(kyvals/ktheta,3)==np.round((kyprime-shell1momentay[k])/ktheta,3))
						indexfind=np.intersect1d(kxfind,kyfind)
						if np.size(indexfind)>0:
							indexnew=indexfind[0]
							kxprime=kxvals[indexnew]+shell1momentax[k]
							kyprime=kyvals[indexnew]+shell1momentay[k]
							Pmatch=np.conjugate(HtotalEigVecsShell1[k][indexnew]).dot(Proj[indexnew]).dot(np.transpose(HtotalEigVecsShell1[k][indexnew]))
					for k in range(np.size(shell2momentax)):
						kxfind=np.where(np.round(kxvals/ktheta,3)==np.round((kxprime-shell2momentax[k])/ktheta,3))
						kyfind=np.where(np.round(kyvals/ktheta,3)==np.round((kyprime-shell2momentay[k])/ktheta,3))
						indexfind=np.intersect1d(kxfind,kyfind)
						if np.size(indexfind)>0:
							indexnew=indexfind[0]
							kxprime=kxvals[indexnew]+shell2momentax[k]
							kyprime=kyvals[indexnew]+shell2momentay[k]
							Pmatch=np.conjugate(HtotalEigVecsShell2[k][indexnew]).dot(Proj[indexnew]).dot(np.transpose(HtotalEigVecsShell2[k][indexnew]))
				else:
					indexnew=indexfind[0]
					kxprime=kxvals[indexnew]
					kyprime=kyvals[indexnew]
					Pmatch=np.conjugate(HtotalEigVecs[indexnew]).dot(Proj[indexnew]).dot(np.transpose(HtotalEigVecs[indexnew]))

				base[i][j]=np.transpose(U).dot(np.conjugate(C2zmatrix).dot(Pmatch).dot(np.transpose(C2zmatrix))).dot(np.conjugate(U))


			elif (kxvals[j]+shell1momentax[i])/ktheta<0+10**(-3) and -(kxvals[j]+shell1momentax[i])/ktheta*q3y/q3x-(kyvals[j]+shell1momentay[i])/ktheta>10**(-3):

				kxprime,kyprime=C3rotate((kxvals[j]+shell1momentax[i]),(kyvals[j]+shell1momentay[i]))
				kxprime,kyprime=C3rotate(kxprime,kyprime)

				kxfind=np.where(np.round(kxvals/ktheta,3)==np.round((kxprime)/ktheta,3))
				kyfind=np.where(np.round(kyvals/ktheta,3)==np.round((kyprime)/ktheta,3))
				indexfind=np.intersect1d(kxfind,kyfind)

				if np.size(indexfind)==0:
					for k in range(np.size(shell1momentax)):
						kxfind=np.where(np.round(kxvals/ktheta,3)==np.round((kxprime-shell1momentax[k])/ktheta,3))
						kyfind=np.where(np.round(kyvals/ktheta,3)==np.round((kyprime-shell1momentay[k])/ktheta,3))
						indexfind=np.intersect1d(kxfind,kyfind)
						if np.size(indexfind)>0:
							indexnew=indexfind[0]
							kxprime=kxvals[indexnew]+shell1momentax[k]
							kyprime=kyvals[indexnew]+shell1momentay[k]
							Pmatch=np.conjugate(HtotalEigVecsShell1[k][indexnew]).dot(Proj[indexnew]).dot(np.transpose(HtotalEigVecsShell1[k][indexnew]))
					for k in range(np.size(shell2momentax)):
						kxfind=np.where(np.round(kxvals/ktheta,3)==np.round((kxprime-shell2momentax[k])/ktheta,3))
						kyfind=np.where(np.round(kyvals/ktheta,3)==np.round((kyprime-shell2momentay[k])/ktheta,3))
						indexfind=np.intersect1d(kxfind,kyfind)
						if np.size(indexfind)>0:
							indexnew=indexfind[0]
							kxprime=kxvals[indexnew]+shell2momentax[k]
							kyprime=kyvals[indexnew]+shell2momentay[k]
							Pmatch=np.conjugate(HtotalEigVecsShell2[k][indexnew]).dot(Proj[indexnew]).dot(np.transpose(HtotalEigVecsShell2[k][indexnew]))
				else:
					indexnew=indexfind[0]
					kxprime=kxvals[indexnew]
					kyprime=kyvals[indexnew]
					Pmatch=np.conjugate(HtotalEigVecs[indexnew]).dot(Proj[indexnew]).dot(np.transpose(HtotalEigVecs[indexnew]))

				base[i][j]=np.transpose(U).dot(np.conjugate(C3matrix).dot(Pmatch).dot(np.transpose(C3matrix))).dot(np.conjugate(U))

			elif (kxvals[j]+shell1momentax[i])/ktheta<0-10**(-3) and -(kxvals[j]+shell1momentax[i])/ktheta*q3y/q3x-(kyvals[j]+shell1momentay[i])/ktheta<10**(-3) and (kxvals[j]+shell1momentax[i])/ktheta*q3y/q3x-(kyvals[j]+shell1momentay[i])/ktheta>10**(-3):
				kxprime,kyprime=C3rotate((kxvals[j]+shell1momentax[i]),(kyvals[j]+shell1momentay[i]))
				kxprime,kyprime=C2zrotate(kxprime,kyprime)

				kxfind=np.where(np.round(kxvals/ktheta,3)==np.round((kxprime)/ktheta,3))
				kyfind=np.where(np.round(kyvals/ktheta,3)==np.round((kyprime)/ktheta,3))
				indexfind=np.intersect1d(kxfind,kyfind)

				if np.size(indexfind)==0:
					for k in range(np.size(shell1momentax)):
						kxfind=np.where(np.round(kxvals/ktheta,3)==np.round((kxprime-shell1momentax[k])/ktheta,3))
						kyfind=np.where(np.round(kyvals/ktheta,3)==np.round((kyprime-shell1momentay[k])/ktheta,3))
						indexfind=np.intersect1d(kxfind,kyfind)
						if np.size(indexfind)>0:
							indexnew=indexfind[0]
							kxprime=kxvals[indexnew]+shell1momentax[k]
							kyprime=kyvals[indexnew]+shell1momentay[k]
							Pmatch=np.conjugate(HtotalEigVecsShell1[k][indexnew]).dot(Proj[indexnew]).dot(np.transpose(HtotalEigVecsShell1[k][indexnew]))
					for k in range(np.size(shell2momentax)):
						kxfind=np.where(np.round(kxvals/ktheta,3)==np.round((kxprime-shell2momentax[k])/ktheta,3))
						kyfind=np.where(np.round(kyvals/ktheta,3)==np.round((kyprime-shell2momentay[k])/ktheta,3))
						indexfind=np.intersect1d(kxfind,kyfind)
						if np.size(indexfind)>0:
							indexnew=indexfind[0]
							kxprime=kxvals[indexnew]+shell2momentax[k]
							kyprime=kyvals[indexnew]+shell2momentay[k]
							Pmatch=np.conjugate(HtotalEigVecsShell2[k][indexnew]).dot(Proj[indexnew]).dot(np.transpose(HtotalEigVecsShell2[k][indexnew]))
				else:
					indexnew=indexfind[0]
					kxprime=kxvals[indexnew]
					kyprime=kyvals[indexnew]
					Pmatch=np.conjugate(HtotalEigVecs[indexnew]).dot(Proj[indexnew]).dot(np.transpose(HtotalEigVecs[indexnew]))
				base[i][j]=np.transpose(U).dot(np.conjugate(C3C3C2product)).dot(Pmatch).dot(np.transpose(C3C3C2product)).dot(np.conjugate(U))
			else:
				print('something bad happened')
	return base
def MakeShell1PC2z(Proj):
	base=np.full((np.size(shell1momentax),npts,nactive*4,nactive*4),0.+0.j)



	for i in range(np.size(shell1momentax)):

		for j in range(npts):
			if not(np.abs(kxvals[j]+b1x)<10**(-3)*ktheta or np.abs(kyvals[j]-kxvals[j]*q2y/q2x-ktheta)<10**(-3) or np.abs(kyvals[j]+kxvals[j]*q2y/q2x-ktheta)<10**(-3)):
				continue
			U=HtotalEigVecsShell1[i][j]


			

			if (kxvals[j]+shell1momentax[i])/ktheta<0-10**(-3) and (kxvals[j]+shell1momentax[i])/ktheta*q3y/q3x-(kyvals[j]+shell1momentay[i])/ktheta<10**(-3):



				base[i][j]=np.transpose(U).dot(np.conjugate(HtotalEigVecsShell1[i][j]).dot(Proj[j]).dot(np.transpose(HtotalEigVecsShell1[i][j]))).dot(np.conjugate(U))

				


			elif (kxvals[j]+shell1momentax[i])/ktheta>0-10**(-3) and -(kxvals[j]+shell1momentax[i])/ktheta*q3y/q3x-(kyvals[j]+shell1momentay[i])/ktheta<-10**(-3):



				kxprime,kyprime=C3rotate((kxvals[j]+shell1momentax[i]),(kyvals[j]+shell1momentay[i]))
				kxprime,kyprime=C3rotate(kxprime,kyprime)
				kxprime,kyprime=C2zrotate(kxprime,kyprime)

				kxfind=np.where(np.round(kxvals/ktheta,3)==np.round((kxprime)/ktheta,3))
				kyfind=np.where(np.round(kyvals/ktheta,3)==np.round((kyprime)/ktheta,3))
				indexfind=np.intersect1d(kxfind,kyfind)

				if np.size(indexfind)==0:
					for k in range(np.size(shell1momentax)):
						kxfind=np.where(np.round(kxvals/ktheta,3)==np.round((kxprime-shell1momentax[k])/ktheta,3))
						kyfind=np.where(np.round(kyvals/ktheta,3)==np.round((kyprime-shell1momentay[k])/ktheta,3))
						indexfind=np.intersect1d(kxfind,kyfind)
						if np.size(indexfind)>0:
							indexnew=indexfind[0]
							kxprime=kxvals[indexnew]+shell1momentax[k]
							kyprime=kyvals[indexnew]+shell1momentay[k]
							Pmatch=np.conjugate(HtotalEigVecsShell1[k][indexnew]).dot(Proj[indexnew]).dot(np.transpose(HtotalEigVecsShell1[k][indexnew]))
					for k in range(np.size(shell2momentax)):
						kxfind=np.where(np.round(kxvals/ktheta,3)==np.round((kxprime-shell2momentax[k])/ktheta,3))
						kyfind=np.where(np.round(kyvals/ktheta,3)==np.round((kyprime-shell2momentay[k])/ktheta,3))
						indexfind=np.intersect1d(kxfind,kyfind)
						if np.size(indexfind)>0:
							indexnew=indexfind[0]
							kxprime=kxvals[indexnew]+shell2momentax[k]
							kyprime=kyvals[indexnew]+shell2momentay[k]
							Pmatch=np.conjugate(HtotalEigVecsShell2[k][indexnew]).dot(Proj[indexnew]).dot(np.transpose(HtotalEigVecsShell2[k][indexnew]))
				else:
					indexnew=indexfind[0]
					kxprime=kxvals[indexnew]
					kyprime=kyvals[indexnew]
					Pmatch=np.conjugate(HtotalEigVecs[indexnew]).dot(Proj[indexnew]).dot(np.transpose(HtotalEigVecs[indexnew]))

				base[i][j]=np.transpose(U).dot(np.conjugate(C3T0product).dot(np.conjugate(Pmatch)).dot(np.transpose(C3T0product))).dot(np.conjugate(U))



			elif (kxvals[j]+shell1momentax[i])/ktheta>0+10**(-3) and -(kxvals[j]+shell1momentax[i])/ktheta*q3y/q3x-(kyvals[j]+shell1momentay[i])/ktheta>-10**(-3) and (kxvals[j]+shell1momentax[i])/ktheta*q3y/q3x-(kyvals[j]+shell1momentay[i])/ktheta<-10**(-3):

				kxprime,kyprime=C3rotate((kxvals[j]+shell1momentax[i]),(kyvals[j]+shell1momentay[i]))

				kxfind=np.where(np.round(kxvals/ktheta,3)==np.round((kxprime)/ktheta,3))
				kyfind=np.where(np.round(kyvals/ktheta,3)==np.round((kyprime)/ktheta,3))
				indexfind=np.intersect1d(kxfind,kyfind)

				if np.size(indexfind)==0:
					for k in range(np.size(shell1momentax)):
						kxfind=np.where(np.round(kxvals/ktheta,3)==np.round((kxprime-shell1momentax[k])/ktheta,3))
						kyfind=np.where(np.round(kyvals/ktheta,3)==np.round((kyprime-shell1momentay[k])/ktheta,3))
						indexfind=np.intersect1d(kxfind,kyfind)
						if np.size(indexfind)>0:
							indexnew=indexfind[0]
							kxprime=kxvals[indexnew]+shell1momentax[k]
							kyprime=kyvals[indexnew]+shell1momentay[k]
							Pmatch=np.conjugate(HtotalEigVecsShell1[k][indexnew]).dot(Proj[indexnew]).dot(np.transpose(HtotalEigVecsShell1[k][indexnew]))
					for k in range(np.size(shell2momentax)):
						kxfind=np.where(np.round(kxvals/ktheta,3)==np.round((kxprime-shell2momentax[k])/ktheta,3))
						kyfind=np.where(np.round(kyvals/ktheta,3)==np.round((kyprime-shell2momentay[k])/ktheta,3))
						indexfind=np.intersect1d(kxfind,kyfind)
						if np.size(indexfind)>0:
							indexnew=indexfind[0]
							kxprime=kxvals[indexnew]+shell2momentax[k]
							kyprime=kyvals[indexnew]+shell2momentay[k]
							Pmatch=np.conjugate(HtotalEigVecsShell2[k][indexnew]).dot(Proj[indexnew]).dot(np.transpose(HtotalEigVecsShell2[k][indexnew]))
				else:
					indexnew=indexfind[0]
					kxprime=kxvals[indexnew]
					kyprime=kyvals[indexnew]
					Pmatch=np.conjugate(HtotalEigVecs[indexnew]).dot(Proj[indexnew]).dot(np.transpose(HtotalEigVecs[indexnew]))

				base[i][j]=np.transpose(U).dot(np.conjugate(C3C3product).dot(Pmatch).dot(np.transpose(C3C3product))).dot(np.conjugate(U))

			elif (kxvals[j]+shell1momentax[i])/ktheta>0+10**(-3) and (kxvals[j]+shell1momentax[i])/ktheta*q3y/q3x-(kyvals[j]+shell1momentay[i])/ktheta>-10**(-3):

				kxprime,kyprime=C2zrotate((kxvals[j]+shell1momentax[i]),(kyvals[j]+shell1momentay[i]))

				kxfind=np.where(np.round(kxvals/ktheta,3)==np.round((kxprime)/ktheta,3))
				kyfind=np.where(np.round(kyvals/ktheta,3)==np.round((kyprime)/ktheta,3))
				indexfind=np.intersect1d(kxfind,kyfind)

				if np.size(indexfind)==0:
					for k in range(np.size(shell1momentax)):
						kxfind=np.where(np.round(kxvals/ktheta,3)==np.round((kxprime-shell1momentax[k])/ktheta,3))
						kyfind=np.where(np.round(kyvals/ktheta,3)==np.round((kyprime-shell1momentay[k])/ktheta,3))
						indexfind=np.intersect1d(kxfind,kyfind)
						if np.size(indexfind)>0:
							indexnew=indexfind[0]
							kxprime=kxvals[indexnew]+shell1momentax[k]
							kyprime=kyvals[indexnew]+shell1momentay[k]
							Pmatch=np.conjugate(HtotalEigVecsShell1[k][indexnew]).dot(Proj[indexnew]).dot(np.transpose(HtotalEigVecsShell1[k][indexnew]))
					for k in range(np.size(shell2momentax)):
						kxfind=np.where(np.round(kxvals/ktheta,3)==np.round((kxprime-shell2momentax[k])/ktheta,3))
						kyfind=np.where(np.round(kyvals/ktheta,3)==np.round((kyprime-shell2momentay[k])/ktheta,3))
						indexfind=np.intersect1d(kxfind,kyfind)
						if np.size(indexfind)>0:
							indexnew=indexfind[0]
							kxprime=kxvals[indexnew]+shell2momentax[k]
							kyprime=kyvals[indexnew]+shell2momentay[k]
							Pmatch=np.conjugate(HtotalEigVecsShell2[k][indexnew]).dot(Proj[indexnew]).dot(np.transpose(HtotalEigVecsShell2[k][indexnew]))
				else:
					indexnew=indexfind[0]
					kxprime=kxvals[indexnew]
					kyprime=kyvals[indexnew]
					Pmatch=np.conjugate(HtotalEigVecs[indexnew]).dot(Proj[indexnew]).dot(np.transpose(HtotalEigVecs[indexnew]))

				base[i][j]=np.transpose(U).dot(np.conjugate(T0matrix).dot(np.conjugate(Pmatch)).dot(np.transpose(T0matrix))).dot(np.conjugate(U))


			elif (kxvals[j]+shell1momentax[i])/ktheta<0+10**(-3) and -(kxvals[j]+shell1momentax[i])/ktheta*q3y/q3x-(kyvals[j]+shell1momentay[i])/ktheta>10**(-3):

				kxprime,kyprime=C3rotate((kxvals[j]+shell1momentax[i]),(kyvals[j]+shell1momentay[i]))
				kxprime,kyprime=C3rotate(kxprime,kyprime)

				kxfind=np.where(np.round(kxvals/ktheta,3)==np.round((kxprime)/ktheta,3))
				kyfind=np.where(np.round(kyvals/ktheta,3)==np.round((kyprime)/ktheta,3))
				indexfind=np.intersect1d(kxfind,kyfind)

				if np.size(indexfind)==0:
					for k in range(np.size(shell1momentax)):
						kxfind=np.where(np.round(kxvals/ktheta,3)==np.round((kxprime-shell1momentax[k])/ktheta,3))
						kyfind=np.where(np.round(kyvals/ktheta,3)==np.round((kyprime-shell1momentay[k])/ktheta,3))
						indexfind=np.intersect1d(kxfind,kyfind)
						if np.size(indexfind)>0:
							indexnew=indexfind[0]
							kxprime=kxvals[indexnew]+shell1momentax[k]
							kyprime=kyvals[indexnew]+shell1momentay[k]
							Pmatch=np.conjugate(HtotalEigVecsShell1[k][indexnew]).dot(Proj[indexnew]).dot(np.transpose(HtotalEigVecsShell1[k][indexnew]))
					for k in range(np.size(shell2momentax)):
						kxfind=np.where(np.round(kxvals/ktheta,3)==np.round((kxprime-shell2momentax[k])/ktheta,3))
						kyfind=np.where(np.round(kyvals/ktheta,3)==np.round((kyprime-shell2momentay[k])/ktheta,3))
						indexfind=np.intersect1d(kxfind,kyfind)
						if np.size(indexfind)>0:
							indexnew=indexfind[0]
							kxprime=kxvals[indexnew]+shell2momentax[k]
							kyprime=kyvals[indexnew]+shell2momentay[k]
							Pmatch=np.conjugate(HtotalEigVecsShell2[k][indexnew]).dot(Proj[indexnew]).dot(np.transpose(HtotalEigVecsShell2[k][indexnew]))
				else:
					indexnew=indexfind[0]
					kxprime=kxvals[indexnew]
					kyprime=kyvals[indexnew]
					Pmatch=np.conjugate(HtotalEigVecs[indexnew]).dot(Proj[indexnew]).dot(np.transpose(HtotalEigVecs[indexnew]))

				base[i][j]=np.transpose(U).dot(np.conjugate(C3matrix).dot(Pmatch).dot(np.transpose(C3matrix))).dot(np.conjugate(U))

			elif (kxvals[j]+shell1momentax[i])/ktheta<0-10**(-3) and -(kxvals[j]+shell1momentax[i])/ktheta*q3y/q3x-(kyvals[j]+shell1momentay[i])/ktheta<10**(-3) and (kxvals[j]+shell1momentax[i])/ktheta*q3y/q3x-(kyvals[j]+shell1momentay[i])/ktheta>10**(-3):
				kxprime,kyprime=C3rotate((kxvals[j]+shell1momentax[i]),(kyvals[j]+shell1momentay[i]))
				kxprime,kyprime=C2zrotate(kxprime,kyprime)

				kxfind=np.where(np.round(kxvals/ktheta,3)==np.round((kxprime)/ktheta,3))
				kyfind=np.where(np.round(kyvals/ktheta,3)==np.round((kyprime)/ktheta,3))
				indexfind=np.intersect1d(kxfind,kyfind)

				if np.size(indexfind)==0:
					for k in range(np.size(shell1momentax)):
						kxfind=np.where(np.round(kxvals/ktheta,3)==np.round((kxprime-shell1momentax[k])/ktheta,3))
						kyfind=np.where(np.round(kyvals/ktheta,3)==np.round((kyprime-shell1momentay[k])/ktheta,3))
						indexfind=np.intersect1d(kxfind,kyfind)
						if np.size(indexfind)>0:
							indexnew=indexfind[0]
							kxprime=kxvals[indexnew]+shell1momentax[k]
							kyprime=kyvals[indexnew]+shell1momentay[k]
							Pmatch=np.conjugate(HtotalEigVecsShell1[k][indexnew]).dot(Proj[indexnew]).dot(np.transpose(HtotalEigVecsShell1[k][indexnew]))
					for k in range(np.size(shell2momentax)):
						kxfind=np.where(np.round(kxvals/ktheta,3)==np.round((kxprime-shell2momentax[k])/ktheta,3))
						kyfind=np.where(np.round(kyvals/ktheta,3)==np.round((kyprime-shell2momentay[k])/ktheta,3))
						indexfind=np.intersect1d(kxfind,kyfind)
						if np.size(indexfind)>0:
							indexnew=indexfind[0]
							kxprime=kxvals[indexnew]+shell2momentax[k]
							kyprime=kyvals[indexnew]+shell2momentay[k]
							Pmatch=np.conjugate(HtotalEigVecsShell2[k][indexnew]).dot(Proj[indexnew]).dot(np.transpose(HtotalEigVecsShell2[k][indexnew]))
				else:
					indexnew=indexfind[0]
					kxprime=kxvals[indexnew]
					kyprime=kyvals[indexnew]
					Pmatch=np.conjugate(HtotalEigVecs[indexnew]).dot(Proj[indexnew]).dot(np.transpose(HtotalEigVecs[indexnew]))
				# print(i,j)

				base[i][j]=np.transpose(U).dot(np.conjugate(C3C3T0product).dot(np.conjugate(Pmatch)).dot(np.transpose(C3C3T0product))).dot(np.conjugate(U))
			else:
				print('something bad happened')


	return base
def MakeShell2P(Proj):


	base=np.full((np.size(shell2momentax),npts,nactive*4,nactive*4),0.+0.j)

	for i in range(np.size(shell2momentax)):
		time0=time.time()

		for j in range(npts):
			if not(np.abs(kxvals[j]+b1x)<10**(-3)*ktheta or np.abs(kyvals[j]-kxvals[j]*q2y/q2x-ktheta)<10**(-3) or np.abs(kyvals[j]+kxvals[j]*q2y/q2x-ktheta)<10**(-3)):
				continue

			if np.sqrt((shell2momentax[i]+kxvals[j])**2+(shell2momentay[i]+kyvals[j])**2)/ktheta >np.sqrt(7+3*np.sqrt(3))-.01 and (np.abs(kxvals[j]+b1x)<10**(-3)*ktheta or np.abs(kyvals[j]-kxvals[j]*q2y/q2x-ktheta)<10**(-3) or np.abs(kyvals[j]+kxvals[j]*q2y/q2x-ktheta)<10**(-3)):
					continue
			U=HtotalEigVecsShell2[i][j]


			

			if (kxvals[j]+shell2momentax[i])/ktheta<0-10**(-3) and (kxvals[j]+shell2momentax[i])/ktheta*q3y/q3x-(kyvals[j]+shell2momentay[i])/ktheta<10**(-3):



				base[i][j]=np.transpose(U).dot(np.conjugate(HtotalEigVecsShell2[i][j]).dot(Proj[j]).dot(np.transpose(HtotalEigVecsShell2[i][j]))).dot(np.conjugate(U))

				


			elif (kxvals[j]+shell2momentax[i])/ktheta>0-10**(-3) and -(kxvals[j]+shell2momentax[i])/ktheta*q3y/q3x-(kyvals[j]+shell2momentay[i])/ktheta<-10**(-3):



				kxprime,kyprime=C3rotate((kxvals[j]+shell2momentax[i]),(kyvals[j]+shell2momentay[i]))
				kxprime,kyprime=C3rotate(kxprime,kyprime)
				kxprime,kyprime=C2zrotate(kxprime,kyprime)

				kxfind=np.where(np.round(kxvals/ktheta,3)==np.round((kxprime)/ktheta,3))
				kyfind=np.where(np.round(kyvals/ktheta,3)==np.round((kyprime)/ktheta,3))
				indexfind=np.intersect1d(kxfind,kyfind)

				if np.size(indexfind)==0:
					for k in range(np.size(shell1momentax)):
						kxfind=np.where(np.round(kxvals/ktheta,3)==np.round((kxprime-shell1momentax[k])/ktheta,3))
						kyfind=np.where(np.round(kyvals/ktheta,3)==np.round((kyprime-shell1momentay[k])/ktheta,3))
						indexfind=np.intersect1d(kxfind,kyfind)
						if np.size(indexfind)>0:
							indexnew=indexfind[0]
							kxprime=kxvals[indexnew]+shell1momentax[k]
							kyprime=kyvals[indexnew]+shell1momentay[k]
							Pmatch=np.conjugate(HtotalEigVecsShell1[k][indexnew]).dot(Proj[indexnew]).dot(np.transpose(HtotalEigVecsShell1[k][indexnew]))
					for k in range(np.size(shell2momentax)):
						kxfind=np.where(np.round(kxvals/ktheta,3)==np.round((kxprime-shell2momentax[k])/ktheta,3))
						kyfind=np.where(np.round(kyvals/ktheta,3)==np.round((kyprime-shell2momentay[k])/ktheta,3))
						indexfind=np.intersect1d(kxfind,kyfind)
						if np.size(indexfind)>0:
							indexnew=indexfind[0]
							kxprime=kxvals[indexnew]+shell2momentax[k]
							kyprime=kyvals[indexnew]+shell2momentay[k]
							Pmatch=np.conjugate(HtotalEigVecsShell2[k][indexnew]).dot(Proj[indexnew]).dot(np.transpose(HtotalEigVecsShell2[k][indexnew]))
					for k in range(np.size(shell3momentax)):
						kxfind=np.where(np.round(kxvals/ktheta,3)==np.round((kxprime-shell3momentax[k])/ktheta,3))
						kyfind=np.where(np.round(kyvals/ktheta,3)==np.round((kyprime-shell3momentay[k])/ktheta,3))
						indexfind=np.intersect1d(kxfind,kyfind)
						if np.size(indexfind)>0:
							indexnew=indexfind[0]
							kxprime=kxvals[indexnew]+shell3momentax[k]
							kyprime=kyvals[indexnew]+shell3momentay[k]
							Pmatch=np.conjugate(HtotalEigVecs[indexnew]).dot(Proj[indexnew]).dot(np.transpose(HtotalEigVecs[indexnew]))

				else:
					indexnew=indexfind[0]
					kxprime=kxvals[indexnew]
					kyprime=kyvals[indexnew]
					Pmatch=np.conjugate(HtotalEigVecs[indexnew]).dot(Proj[indexnew]).dot(np.transpose(HtotalEigVecs[indexnew]))
				time1=time.time()

				base[i][j]=np.transpose(U).dot(np.conjugate(C3C2product).dot(Pmatch).dot(np.transpose(C3C2product))).dot(np.conjugate(U))




			elif (kxvals[j]+shell2momentax[i])/ktheta>0+10**(-3) and -(kxvals[j]+shell2momentax[i])/ktheta*q3y/q3x-(kyvals[j]+shell2momentay[i])/ktheta>-10**(-3) and (kxvals[j]+shell2momentax[i])/ktheta*q3y/q3x-(kyvals[j]+shell2momentay[i])/ktheta<-10**(-3):

				kxprime,kyprime=C3rotate((kxvals[j]+shell2momentax[i]),(kyvals[j]+shell2momentay[i]))

				kxfind=np.where(np.round(kxvals/ktheta,3)==np.round((kxprime)/ktheta,3))
				kyfind=np.where(np.round(kyvals/ktheta,3)==np.round((kyprime)/ktheta,3))
				indexfind=np.intersect1d(kxfind,kyfind)

				if np.size(indexfind)==0:
					for k in range(np.size(shell1momentax)):
						kxfind=np.where(np.round(kxvals/ktheta,3)==np.round((kxprime-shell1momentax[k])/ktheta,3))
						kyfind=np.where(np.round(kyvals/ktheta,3)==np.round((kyprime-shell1momentay[k])/ktheta,3))
						indexfind=np.intersect1d(kxfind,kyfind)
						if np.size(indexfind)>0:
							indexnew=indexfind[0]
							kxprime=kxvals[indexnew]+shell1momentax[k]
							kyprime=kyvals[indexnew]+shell1momentay[k]
							Pmatch=np.conjugate(HtotalEigVecsShell1[k][indexnew]).dot(Proj[indexnew]).dot(np.transpose(HtotalEigVecsShell1[k][indexnew]))
					for k in range(np.size(shell2momentax)):
						kxfind=np.where(np.round(kxvals/ktheta,3)==np.round((kxprime-shell2momentax[k])/ktheta,3))
						kyfind=np.where(np.round(kyvals/ktheta,3)==np.round((kyprime-shell2momentay[k])/ktheta,3))
						indexfind=np.intersect1d(kxfind,kyfind)
						if np.size(indexfind)>0:
							indexnew=indexfind[0]
							kxprime=kxvals[indexnew]+shell2momentax[k]
							kyprime=kyvals[indexnew]+shell2momentay[k]
							Pmatch=np.conjugate(HtotalEigVecsShell2[k][indexnew]).dot(Proj[indexnew]).dot(np.transpose(HtotalEigVecsShell2[k][indexnew]))

					for k in range(np.size(shell3momentax)):
						kxfind=np.where(np.round(kxvals/ktheta,3)==np.round((kxprime-shell3momentax[k])/ktheta,3))
						kyfind=np.where(np.round(kyvals/ktheta,3)==np.round((kyprime-shell3momentay[k])/ktheta,3))
						indexfind=np.intersect1d(kxfind,kyfind)
						if np.size(indexfind)>0:
							indexnew=indexfind[0]
							kxprime=kxvals[indexnew]+shell3momentax[k]
							kyprime=kyvals[indexnew]+shell3momentay[k]
							Pmatch=np.conjugate(HtotalEigVecs[indexnew]).dot(Proj[indexnew]).dot(np.transpose(HtotalEigVecs[indexnew]))
				else:
					indexnew=indexfind[0]
					kxprime=kxvals[indexnew]
					kyprime=kyvals[indexnew]
					Pmatch=np.conjugate(HtotalEigVecs[indexnew]).dot(Proj[indexnew]).dot(np.transpose(HtotalEigVecs[indexnew]))

				base[i][j]=np.transpose(U).dot(np.conjugate(C3C3product).dot(Pmatch).dot(np.transpose(C3C3product))).dot(np.conjugate(U))

			elif (kxvals[j]+shell2momentax[i])/ktheta>0+10**(-3) and (kxvals[j]+shell2momentax[i])/ktheta*q3y/q3x-(kyvals[j]+shell2momentay[i])/ktheta>-10**(-3):

				kxprime,kyprime=C2zrotate((kxvals[j]+shell2momentax[i]),(kyvals[j]+shell2momentay[i]))

				kxfind=np.where(np.round(kxvals/ktheta,3)==np.round((kxprime)/ktheta,3))
				kyfind=np.where(np.round(kyvals/ktheta,3)==np.round((kyprime)/ktheta,3))
				indexfind=np.intersect1d(kxfind,kyfind)

				if np.size(indexfind)==0:
					for k in range(np.size(shell1momentax)):
						kxfind=np.where(np.round(kxvals/ktheta,3)==np.round((kxprime-shell1momentax[k])/ktheta,3))
						kyfind=np.where(np.round(kyvals/ktheta,3)==np.round((kyprime-shell1momentay[k])/ktheta,3))
						indexfind=np.intersect1d(kxfind,kyfind)
						if np.size(indexfind)>0:
							indexnew=indexfind[0]
							kxprime=kxvals[indexnew]+shell1momentax[k]
							kyprime=kyvals[indexnew]+shell1momentay[k]
							Pmatch=np.conjugate(HtotalEigVecsShell1[k][indexnew]).dot(Proj[indexnew]).dot(np.transpose(HtotalEigVecsShell1[k][indexnew]))
					for k in range(np.size(shell2momentax)):
						kxfind=np.where(np.round(kxvals/ktheta,3)==np.round((kxprime-shell2momentax[k])/ktheta,3))
						kyfind=np.where(np.round(kyvals/ktheta,3)==np.round((kyprime-shell2momentay[k])/ktheta,3))
						indexfind=np.intersect1d(kxfind,kyfind)
						if np.size(indexfind)>0:
							indexnew=indexfind[0]
							kxprime=kxvals[indexnew]+shell2momentax[k]
							kyprime=kyvals[indexnew]+shell2momentay[k]
							Pmatch=np.conjugate(HtotalEigVecsShell2[k][indexnew]).dot(Proj[indexnew]).dot(np.transpose(HtotalEigVecsShell2[k][indexnew]))
					for k in range(np.size(shell3momentax)):
						kxfind=np.where(np.round(kxvals/ktheta,3)==np.round((kxprime-shell3momentax[k])/ktheta,3))
						kyfind=np.where(np.round(kyvals/ktheta,3)==np.round((kyprime-shell3momentay[k])/ktheta,3))
						indexfind=np.intersect1d(kxfind,kyfind)
						if np.size(indexfind)>0:
							indexnew=indexfind[0]
							kxprime=kxvals[indexnew]+shell3momentax[k]
							kyprime=kyvals[indexnew]+shell3momentay[k]
							Pmatch=np.conjugate(HtotalEigVecs[indexnew]).dot(Proj[indexnew]).dot(np.transpose(HtotalEigVecs[indexnew]))
				else:
					indexnew=indexfind[0]
					kxprime=kxvals[indexnew]
					kyprime=kyvals[indexnew]
					Pmatch=np.conjugate(HtotalEigVecs[indexnew]).dot(Proj[indexnew]).dot(np.transpose(HtotalEigVecs[indexnew]))

				base[i][j]=np.transpose(U).dot(np.conjugate(C2zmatrix).dot(Pmatch).dot(np.transpose(C2zmatrix))).dot(np.conjugate(U))


			elif (kxvals[j]+shell2momentax[i])/ktheta<0+10**(-3) and -(kxvals[j]+shell2momentax[i])/ktheta*q3y/q3x-(kyvals[j]+shell2momentay[i])/ktheta>10**(-3):

				kxprime,kyprime=C3rotate((kxvals[j]+shell2momentax[i]),(kyvals[j]+shell2momentay[i]))
				kxprime,kyprime=C3rotate(kxprime,kyprime)

				kxfind=np.where(np.round(kxvals/ktheta,3)==np.round((kxprime)/ktheta,3))
				kyfind=np.where(np.round(kyvals/ktheta,3)==np.round((kyprime)/ktheta,3))
				indexfind=np.intersect1d(kxfind,kyfind)

				if np.size(indexfind)==0:
					for k in range(np.size(shell1momentax)):
						kxfind=np.where(np.round(kxvals/ktheta,3)==np.round((kxprime-shell1momentax[k])/ktheta,3))
						kyfind=np.where(np.round(kyvals/ktheta,3)==np.round((kyprime-shell1momentay[k])/ktheta,3))
						indexfind=np.intersect1d(kxfind,kyfind)
						if np.size(indexfind)>0:
							indexnew=indexfind[0]
							kxprime=kxvals[indexnew]+shell1momentax[k]
							kyprime=kyvals[indexnew]+shell1momentay[k]
							Pmatch=np.conjugate(HtotalEigVecsShell1[k][indexnew]).dot(Proj[indexnew]).dot(np.transpose(HtotalEigVecsShell1[k][indexnew]))
					for k in range(np.size(shell2momentax)):
						kxfind=np.where(np.round(kxvals/ktheta,3)==np.round((kxprime-shell2momentax[k])/ktheta,3))
						kyfind=np.where(np.round(kyvals/ktheta,3)==np.round((kyprime-shell2momentay[k])/ktheta,3))
						indexfind=np.intersect1d(kxfind,kyfind)
						if np.size(indexfind)>0:
							indexnew=indexfind[0]
							kxprime=kxvals[indexnew]+shell2momentax[k]
							kyprime=kyvals[indexnew]+shell2momentay[k]
							Pmatch=np.conjugate(HtotalEigVecsShell2[k][indexnew]).dot(Proj[indexnew]).dot(np.transpose(HtotalEigVecsShell2[k][indexnew]))
					for k in range(np.size(shell3momentax)):
						kxfind=np.where(np.round(kxvals/ktheta,3)==np.round((kxprime-shell3momentax[k])/ktheta,3))
						kyfind=np.where(np.round(kyvals/ktheta,3)==np.round((kyprime-shell3momentay[k])/ktheta,3))
						indexfind=np.intersect1d(kxfind,kyfind)
						if np.size(indexfind)>0:
							indexnew=indexfind[0]
							kxprime=kxvals[indexnew]+shell3momentax[k]
							kyprime=kyvals[indexnew]+shell3momentay[k]
							Pmatch=np.conjugate(HtotalEigVecs[indexnew]).dot(Proj[indexnew]).dot(np.transpose(HtotalEigVecs[indexnew]))
				else:
					indexnew=indexfind[0]
					kxprime=kxvals[indexnew]
					kyprime=kyvals[indexnew]
					Pmatch=np.conjugate(HtotalEigVecs[indexnew]).dot(Proj[indexnew]).dot(np.transpose(HtotalEigVecs[indexnew]))

				base[i][j]=np.transpose(U).dot(np.conjugate(C3matrix).dot(Pmatch).dot(np.transpose(C3matrix))).dot(np.conjugate(U))

			elif (kxvals[j]+shell2momentax[i])/ktheta<0-10**(-3) and -(kxvals[j]+shell2momentax[i])/ktheta*q3y/q3x-(kyvals[j]+shell2momentay[i])/ktheta<10**(-3) and (kxvals[j]+shell2momentax[i])/ktheta*q3y/q3x-(kyvals[j]+shell2momentay[i])/ktheta>10**(-3):
				kxprime,kyprime=C3rotate((kxvals[j]+shell2momentax[i]),(kyvals[j]+shell2momentay[i]))
				kxprime,kyprime=C2zrotate(kxprime,kyprime)

				kxfind=np.where(np.round(kxvals/ktheta,3)==np.round((kxprime)/ktheta,3))
				kyfind=np.where(np.round(kyvals/ktheta,3)==np.round((kyprime)/ktheta,3))
				indexfind=np.intersect1d(kxfind,kyfind)

				if np.size(indexfind)==0:
					for k in range(np.size(shell1momentax)):
						kxfind=np.where(np.round(kxvals/ktheta,3)==np.round((kxprime-shell1momentax[k])/ktheta,3))
						kyfind=np.where(np.round(kyvals/ktheta,3)==np.round((kyprime-shell1momentay[k])/ktheta,3))
						indexfind=np.intersect1d(kxfind,kyfind)
						if np.size(indexfind)>0:
							indexnew=indexfind[0]
							kxprime=kxvals[indexnew]+shell1momentax[k]
							kyprime=kyvals[indexnew]+shell1momentay[k]
							Pmatch=np.conjugate(HtotalEigVecsShell1[k][indexnew]).dot(Proj[indexnew]).dot(np.transpose(HtotalEigVecsShell1[k][indexnew]))
					for k in range(np.size(shell2momentax)):
						kxfind=np.where(np.round(kxvals/ktheta,3)==np.round((kxprime-shell2momentax[k])/ktheta,3))
						kyfind=np.where(np.round(kyvals/ktheta,3)==np.round((kyprime-shell2momentay[k])/ktheta,3))
						indexfind=np.intersect1d(kxfind,kyfind)
						if np.size(indexfind)>0:
							indexnew=indexfind[0]
							kxprime=kxvals[indexnew]+shell2momentax[k]
							kyprime=kyvals[indexnew]+shell2momentay[k]
							Pmatch=np.conjugate(HtotalEigVecsShell2[k][indexnew]).dot(Proj[indexnew]).dot(np.transpose(HtotalEigVecsShell2[k][indexnew]))
					for k in range(np.size(shell3momentax)):
						kxfind=np.where(np.round(kxvals/ktheta,3)==np.round((kxprime-shell3momentax[k])/ktheta,3))
						kyfind=np.where(np.round(kyvals/ktheta,3)==np.round((kyprime-shell3momentay[k])/ktheta,3))
						indexfind=np.intersect1d(kxfind,kyfind)
						if np.size(indexfind)>0:
							indexnew=indexfind[0]
							kxprime=kxvals[indexnew]+shell3momentax[k]
							kyprime=kyvals[indexnew]+shell3momentay[k]
				else:
					indexnew=indexfind[0]
					kxprime=kxvals[indexnew]
					kyprime=kyvals[indexnew]
					Pmatch=np.conjugate(HtotalEigVecs[indexnew]).dot(Proj[indexnew]).dot(np.transpose(HtotalEigVecs[indexnew]))

				base[i][j]=np.transpose(U).dot(np.conjugate(C3C3C2product).dot(Pmatch).dot(np.transpose(C3C3C2product))).dot(np.conjugate(U))
			else:
				print('something bad happened')
			U = None
			del U


	return base
def MakeShell2PC2z(Proj):
	base=np.full((np.size(shell2momentax),npts,nactive*4,nactive*4),0.+0.j)

	for i in range(np.size(shell2momentax)):

		for j in range(npts):
			if not(np.abs(kxvals[j]+b1x)<10**(-3)*ktheta or np.abs(kyvals[j]-kxvals[j]*q2y/q2x-ktheta)<10**(-3) or np.abs(kyvals[j]+kxvals[j]*q2y/q2x-ktheta)<10**(-3)):
				continue
			if np.sqrt((shell2momentax[i]+kxvals[j])**2+(shell2momentay[i]+kyvals[j])**2)/ktheta >np.sqrt(7+3*np.sqrt(3))-.01 and (np.abs(kxvals[j]+b1x)<10**(-3)*ktheta or np.abs(kyvals[j]-kxvals[j]*q2y/q2x-ktheta)<10**(-3) or np.abs(kyvals[j]+kxvals[j]*q2y/q2x-ktheta)<10**(-3)):
					continue
			U=HtotalEigVecsShell2[i][j]


			

			if (kxvals[j]+shell2momentax[i])/ktheta<0-10**(-3) and (kxvals[j]+shell2momentax[i])/ktheta*q3y/q3x-(kyvals[j]+shell2momentay[i])/ktheta<10**(-3):



				base[i][j]=np.transpose(U).dot(np.conjugate(HtotalEigVecsShell2[i][j]).dot(Proj[j]).dot(np.transpose(HtotalEigVecsShell2[i][j]))).dot(np.conjugate(U))

				


			elif (kxvals[j]+shell2momentax[i])/ktheta>0-10**(-3) and -(kxvals[j]+shell2momentax[i])/ktheta*q3y/q3x-(kyvals[j]+shell2momentay[i])/ktheta<-10**(-3):



				kxprime,kyprime=C3rotate((kxvals[j]+shell2momentax[i]),(kyvals[j]+shell2momentay[i]))
				kxprime,kyprime=C3rotate(kxprime,kyprime)
				kxprime,kyprime=C2zrotate(kxprime,kyprime)

				kxfind=np.where(np.round(kxvals/ktheta,3)==np.round((kxprime)/ktheta,3))
				kyfind=np.where(np.round(kyvals/ktheta,3)==np.round((kyprime)/ktheta,3))
				indexfind=np.intersect1d(kxfind,kyfind)

				if np.size(indexfind)==0:
					for k in range(np.size(shell1momentax)):
						kxfind=np.where(np.round(kxvals/ktheta,3)==np.round((kxprime-shell1momentax[k])/ktheta,3))
						kyfind=np.where(np.round(kyvals/ktheta,3)==np.round((kyprime-shell1momentay[k])/ktheta,3))
						indexfind=np.intersect1d(kxfind,kyfind)
						if np.size(indexfind)>0:
							indexnew=indexfind[0]
							kxprime=kxvals[indexnew]+shell1momentax[k]
							kyprime=kyvals[indexnew]+shell1momentay[k]
							Pmatch=np.conjugate(HtotalEigVecsShell1[k][indexnew]).dot(Proj[indexnew]).dot(np.transpose(HtotalEigVecsShell1[k][indexnew]))
					for k in range(np.size(shell2momentax)):
						kxfind=np.where(np.round(kxvals/ktheta,3)==np.round((kxprime-shell2momentax[k])/ktheta,3))
						kyfind=np.where(np.round(kyvals/ktheta,3)==np.round((kyprime-shell2momentay[k])/ktheta,3))
						indexfind=np.intersect1d(kxfind,kyfind)
						if np.size(indexfind)>0:
							indexnew=indexfind[0]
							kxprime=kxvals[indexnew]+shell2momentax[k]
							kyprime=kyvals[indexnew]+shell2momentay[k]
							Pmatch=np.conjugate(HtotalEigVecsShell2[k][indexnew]).dot(Proj[indexnew]).dot(np.transpose(HtotalEigVecsShell2[k][indexnew]))
					for k in range(np.size(shell3momentax)):
						kxfind=np.where(np.round(kxvals/ktheta,3)==np.round((kxprime-shell3momentax[k])/ktheta,3))
						kyfind=np.where(np.round(kyvals/ktheta,3)==np.round((kyprime-shell3momentay[k])/ktheta,3))
						indexfind=np.intersect1d(kxfind,kyfind)
						if np.size(indexfind)>0:
							indexnew=indexfind[0]
							kxprime=kxvals[indexnew]+shell3momentax[k]
							kyprime=kyvals[indexnew]+shell3momentay[k]
							Pmatch=np.conjugate(HtotalEigVecs[indexnew]).dot(Proj[indexnew]).dot(np.transpose(HtotalEigVecs[indexnew]))

				else:
					indexnew=indexfind[0]
					kxprime=kxvals[indexnew]
					kyprime=kyvals[indexnew]
					Pmatch=np.conjugate(HtotalEigVecs[indexnew]).dot(Proj[indexnew]).dot(np.transpose(HtotalEigVecs[indexnew]))

				base[i][j]=np.transpose(U).dot(np.conjugate(C3matrix).dot(np.conjugate(T0matrix)).dot(np.conjugate(Pmatch)).dot(np.transpose(T0matrix)).dot(np.transpose(C3matrix))).dot(np.conjugate(U))



			elif (kxvals[j]+shell2momentax[i])/ktheta>0+10**(-3) and -(kxvals[j]+shell2momentax[i])/ktheta*q3y/q3x-(kyvals[j]+shell2momentay[i])/ktheta>-10**(-3) and (kxvals[j]+shell2momentax[i])/ktheta*q3y/q3x-(kyvals[j]+shell2momentay[i])/ktheta<-10**(-3):

				kxprime,kyprime=C3rotate((kxvals[j]+shell2momentax[i]),(kyvals[j]+shell2momentay[i]))

				kxfind=np.where(np.round(kxvals/ktheta,3)==np.round((kxprime)/ktheta,3))
				kyfind=np.where(np.round(kyvals/ktheta,3)==np.round((kyprime)/ktheta,3))
				indexfind=np.intersect1d(kxfind,kyfind)

				if np.size(indexfind)==0:
					for k in range(np.size(shell1momentax)):
						kxfind=np.where(np.round(kxvals/ktheta,3)==np.round((kxprime-shell1momentax[k])/ktheta,3))
						kyfind=np.where(np.round(kyvals/ktheta,3)==np.round((kyprime-shell1momentay[k])/ktheta,3))
						indexfind=np.intersect1d(kxfind,kyfind)
						if np.size(indexfind)>0:
							indexnew=indexfind[0]
							kxprime=kxvals[indexnew]+shell1momentax[k]
							kyprime=kyvals[indexnew]+shell1momentay[k]
							Pmatch=np.conjugate(HtotalEigVecsShell1[k][indexnew]).dot(Proj[indexnew]).dot(np.transpose(HtotalEigVecsShell1[k][indexnew]))
					for k in range(np.size(shell2momentax)):
						kxfind=np.where(np.round(kxvals/ktheta,3)==np.round((kxprime-shell2momentax[k])/ktheta,3))
						kyfind=np.where(np.round(kyvals/ktheta,3)==np.round((kyprime-shell2momentay[k])/ktheta,3))
						indexfind=np.intersect1d(kxfind,kyfind)
						if np.size(indexfind)>0:
							indexnew=indexfind[0]
							kxprime=kxvals[indexnew]+shell2momentax[k]
							kyprime=kyvals[indexnew]+shell2momentay[k]
							Pmatch=np.conjugate(HtotalEigVecsShell2[k][indexnew]).dot(Proj[indexnew]).dot(np.transpose(HtotalEigVecsShell2[k][indexnew]))

					for k in range(np.size(shell3momentax)):
						kxfind=np.where(np.round(kxvals/ktheta,3)==np.round((kxprime-shell3momentax[k])/ktheta,3))
						kyfind=np.where(np.round(kyvals/ktheta,3)==np.round((kyprime-shell3momentay[k])/ktheta,3))
						indexfind=np.intersect1d(kxfind,kyfind)
						if np.size(indexfind)>0:
							indexnew=indexfind[0]
							kxprime=kxvals[indexnew]+shell3momentax[k]
							kyprime=kyvals[indexnew]+shell3momentay[k]
							Pmatch=np.conjugate(HtotalEigVecs[indexnew]).dot(Proj[indexnew]).dot(np.transpose(HtotalEigVecs[indexnew]))
				else:
					indexnew=indexfind[0]
					kxprime=kxvals[indexnew]
					kyprime=kyvals[indexnew]
					Pmatch=np.conjugate(HtotalEigVecs[indexnew]).dot(Proj[indexnew]).dot(np.transpose(HtotalEigVecs[indexnew]))

				base[i][j]=np.transpose(U).dot(np.conjugate(C3matrix).dot(np.conjugate(C3matrix)).dot(Pmatch).dot(np.transpose(C3matrix)).dot(np.transpose(C3matrix))).dot(np.conjugate(U))

			elif (kxvals[j]+shell2momentax[i])/ktheta>0+10**(-3) and (kxvals[j]+shell2momentax[i])/ktheta*q3y/q3x-(kyvals[j]+shell2momentay[i])/ktheta>-10**(-3):

				kxprime,kyprime=C2zrotate((kxvals[j]+shell2momentax[i]),(kyvals[j]+shell2momentay[i]))

				kxfind=np.where(np.round(kxvals/ktheta,3)==np.round((kxprime)/ktheta,3))
				kyfind=np.where(np.round(kyvals/ktheta,3)==np.round((kyprime)/ktheta,3))
				indexfind=np.intersect1d(kxfind,kyfind)

				if np.size(indexfind)==0:
					for k in range(np.size(shell1momentax)):
						kxfind=np.where(np.round(kxvals/ktheta,3)==np.round((kxprime-shell1momentax[k])/ktheta,3))
						kyfind=np.where(np.round(kyvals/ktheta,3)==np.round((kyprime-shell1momentay[k])/ktheta,3))
						indexfind=np.intersect1d(kxfind,kyfind)
						if np.size(indexfind)>0:
							indexnew=indexfind[0]
							kxprime=kxvals[indexnew]+shell1momentax[k]
							kyprime=kyvals[indexnew]+shell1momentay[k]
							Pmatch=np.conjugate(HtotalEigVecsShell1[k][indexnew]).dot(Proj[indexnew]).dot(np.transpose(HtotalEigVecsShell1[k][indexnew]))
					for k in range(np.size(shell2momentax)):
						kxfind=np.where(np.round(kxvals/ktheta,3)==np.round((kxprime-shell2momentax[k])/ktheta,3))
						kyfind=np.where(np.round(kyvals/ktheta,3)==np.round((kyprime-shell2momentay[k])/ktheta,3))
						indexfind=np.intersect1d(kxfind,kyfind)
						if np.size(indexfind)>0:
							indexnew=indexfind[0]
							kxprime=kxvals[indexnew]+shell2momentax[k]
							kyprime=kyvals[indexnew]+shell2momentay[k]
							Pmatch=np.conjugate(HtotalEigVecsShell2[k][indexnew]).dot(Proj[indexnew]).dot(np.transpose(HtotalEigVecsShell2[k][indexnew]))
					for k in range(np.size(shell3momentax)):
						kxfind=np.where(np.round(kxvals/ktheta,3)==np.round((kxprime-shell3momentax[k])/ktheta,3))
						kyfind=np.where(np.round(kyvals/ktheta,3)==np.round((kyprime-shell3momentay[k])/ktheta,3))
						indexfind=np.intersect1d(kxfind,kyfind)
						if np.size(indexfind)>0:
							indexnew=indexfind[0]
							kxprime=kxvals[indexnew]+shell3momentax[k]
							kyprime=kyvals[indexnew]+shell3momentay[k]
							Pmatch=np.conjugate(HtotalEigVecs[indexnew]).dot(Proj[indexnew]).dot(np.transpose(HtotalEigVecs[indexnew]))
				else:
					indexnew=indexfind[0]
					kxprime=kxvals[indexnew]
					kyprime=kyvals[indexnew]
					Pmatch=np.conjugate(HtotalEigVecs[indexnew]).dot(Proj[indexnew]).dot(np.transpose(HtotalEigVecs[indexnew]))

				base[i][j]=np.transpose(U).dot(np.conjugate(T0matrix).dot(np.conjugate(Pmatch)).dot(np.transpose(T0matrix))).dot(np.conjugate(U))


			elif (kxvals[j]+shell2momentax[i])/ktheta<0+10**(-3) and -(kxvals[j]+shell2momentax[i])/ktheta*q3y/q3x-(kyvals[j]+shell2momentay[i])/ktheta>10**(-3):

				kxprime,kyprime=C3rotate((kxvals[j]+shell2momentax[i]),(kyvals[j]+shell2momentay[i]))
				kxprime,kyprime=C3rotate(kxprime,kyprime)

				kxfind=np.where(np.round(kxvals/ktheta,3)==np.round((kxprime)/ktheta,3))
				kyfind=np.where(np.round(kyvals/ktheta,3)==np.round((kyprime)/ktheta,3))
				indexfind=np.intersect1d(kxfind,kyfind)

				if np.size(indexfind)==0:
					for k in range(np.size(shell1momentax)):
						kxfind=np.where(np.round(kxvals/ktheta,3)==np.round((kxprime-shell1momentax[k])/ktheta,3))
						kyfind=np.where(np.round(kyvals/ktheta,3)==np.round((kyprime-shell1momentay[k])/ktheta,3))
						indexfind=np.intersect1d(kxfind,kyfind)
						if np.size(indexfind)>0:
							indexnew=indexfind[0]
							kxprime=kxvals[indexnew]+shell1momentax[k]
							kyprime=kyvals[indexnew]+shell1momentay[k]
							Pmatch=np.conjugate(HtotalEigVecsShell1[k][indexnew]).dot(Proj[indexnew]).dot(np.transpose(HtotalEigVecsShell1[k][indexnew]))
					for k in range(np.size(shell2momentax)):
						kxfind=np.where(np.round(kxvals/ktheta,3)==np.round((kxprime-shell2momentax[k])/ktheta,3))
						kyfind=np.where(np.round(kyvals/ktheta,3)==np.round((kyprime-shell2momentay[k])/ktheta,3))
						indexfind=np.intersect1d(kxfind,kyfind)
						if np.size(indexfind)>0:
							indexnew=indexfind[0]
							kxprime=kxvals[indexnew]+shell2momentax[k]
							kyprime=kyvals[indexnew]+shell2momentay[k]
							Pmatch=np.conjugate(HtotalEigVecsShell2[k][indexnew]).dot(Proj[indexnew]).dot(np.transpose(HtotalEigVecsShell2[k][indexnew]))
					for k in range(np.size(shell3momentax)):
						kxfind=np.where(np.round(kxvals/ktheta,3)==np.round((kxprime-shell3momentax[k])/ktheta,3))
						kyfind=np.where(np.round(kyvals/ktheta,3)==np.round((kyprime-shell3momentay[k])/ktheta,3))
						indexfind=np.intersect1d(kxfind,kyfind)
						if np.size(indexfind)>0:
							indexnew=indexfind[0]
							kxprime=kxvals[indexnew]+shell3momentax[k]
							kyprime=kyvals[indexnew]+shell3momentay[k]
							Pmatch=np.conjugate(HtotalEigVecs[indexnew]).dot(Proj[indexnew]).dot(np.transpose(HtotalEigVecs[indexnew]))
				else:
					indexnew=indexfind[0]
					kxprime=kxvals[indexnew]
					kyprime=kyvals[indexnew]
					Pmatch=np.conjugate(HtotalEigVecs[indexnew]).dot(Proj[indexnew]).dot(np.transpose(HtotalEigVecs[indexnew]))

				base[i][j]=np.transpose(U).dot(np.conjugate(C3matrix).dot(Pmatch).dot(np.transpose(C3matrix))).dot(np.conjugate(U))

			elif (kxvals[j]+shell2momentax[i])/ktheta<0-10**(-3) and -(kxvals[j]+shell2momentax[i])/ktheta*q3y/q3x-(kyvals[j]+shell2momentay[i])/ktheta<10**(-3) and (kxvals[j]+shell2momentax[i])/ktheta*q3y/q3x-(kyvals[j]+shell2momentay[i])/ktheta>10**(-3):
				kxprime,kyprime=C3rotate((kxvals[j]+shell2momentax[i]),(kyvals[j]+shell2momentay[i]))
				kxprime,kyprime=C2zrotate(kxprime,kyprime)

				kxfind=np.where(np.round(kxvals/ktheta,3)==np.round((kxprime)/ktheta,3))
				kyfind=np.where(np.round(kyvals/ktheta,3)==np.round((kyprime)/ktheta,3))
				indexfind=np.intersect1d(kxfind,kyfind)

				if np.size(indexfind)==0:
					for k in range(np.size(shell1momentax)):
						kxfind=np.where(np.round(kxvals/ktheta,3)==np.round((kxprime-shell1momentax[k])/ktheta,3))
						kyfind=np.where(np.round(kyvals/ktheta,3)==np.round((kyprime-shell1momentay[k])/ktheta,3))
						indexfind=np.intersect1d(kxfind,kyfind)
						if np.size(indexfind)>0:
							indexnew=indexfind[0]
							kxprime=kxvals[indexnew]+shell1momentax[k]
							kyprime=kyvals[indexnew]+shell1momentay[k]
							Pmatch=np.conjugate(HtotalEigVecsShell1[k][indexnew]).dot(Proj[indexnew]).dot(np.transpose(HtotalEigVecsShell1[k][indexnew]))
					for k in range(np.size(shell2momentax)):
						kxfind=np.where(np.round(kxvals/ktheta,3)==np.round((kxprime-shell2momentax[k])/ktheta,3))
						kyfind=np.where(np.round(kyvals/ktheta,3)==np.round((kyprime-shell2momentay[k])/ktheta,3))
						indexfind=np.intersect1d(kxfind,kyfind)
						if np.size(indexfind)>0:
							indexnew=indexfind[0]
							kxprime=kxvals[indexnew]+shell2momentax[k]
							kyprime=kyvals[indexnew]+shell2momentay[k]
							Pmatch=np.conjugate(HtotalEigVecsShell2[k][indexnew]).dot(Proj[indexnew]).dot(np.transpose(HtotalEigVecsShell2[k][indexnew]))
					for k in range(np.size(shell3momentax)):
						kxfind=np.where(np.round(kxvals/ktheta,3)==np.round((kxprime-shell3momentax[k])/ktheta,3))
						kyfind=np.where(np.round(kyvals/ktheta,3)==np.round((kyprime-shell3momentay[k])/ktheta,3))
						indexfind=np.intersect1d(kxfind,kyfind)
						if np.size(indexfind)>0:
							indexnew=indexfind[0]
							kxprime=kxvals[indexnew]+shell3momentax[k]
							kyprime=kyvals[indexnew]+shell3momentay[k]
							Pmatch=np.conjugate(HtotalEigVecs[indexnew]).dot(Proj[indexnew]).dot(np.transpose(HtotalEigVecs[indexnew]))
				else:
					indexnew=indexfind[0]
					kxprime=kxvals[indexnew]
					kyprime=kyvals[indexnew]
					Pmatch=np.conjugate(HtotalEigVecs[indexnew]).dot(Proj[indexnew]).dot(np.transpose(HtotalEigVecs[indexnew]))
				# print(i,j)

				base[i][j]=np.transpose(U).dot(np.conjugate(C3matrix).dot(np.conjugate(C3matrix)).dot(T0matrix).dot(np.conjugate(Pmatch)).dot(T0matrix).dot(np.transpose(C3matrix)).dot(np.transpose(C3matrix))).dot(np.conjugate(U))
			else:
				print('something bad happened')

	return base
# Useful operators
def MZ():
	nbands=np.size(qxvalstot)*2
	U=np.full((nbands,nbands),0.+0.j)
	Uminus=np.full((nbands,nbands),0.+0.j)
	for i in range(np.size(qxvalstot)):
		for j in range(np.size(qxvalstot)):

			if qvalslayer[i]==2 and qvalslayer[j]==2 and qxvalstot[i]==qxvalstot[j] and qyvalstot[i]==qyvalstot[j]:
				U[2*i][2*j]=1
				U[2*i+1][2*j+1]=1


			if qvalslayer[i]==1 and qvalslayer[j]==3 and qxvalstot[i]==qxvalstot[j] and qyvalstot[i]==qyvalstot[j]:
				U[2*i][2*j]=1
				U[2*i+1][2*j+1]=1
			if qvalslayer[i]==3 and qvalslayer[j]==1 and qxvalstot[i]==qxvalstot[j] and qyvalstot[i]==qyvalstot[j]:
				U[2*i][2*j]=1
				U[2*i+1][2*j+1]=1

	for i in range(np.size(qxvalstotminus)):
		for j in range(np.size(qxvalstotminus)):


			if qvalslayerminus[i]==2 and qvalslayerminus[j]==2 and qxvalstotminus[i]==qxvalstotminus[j] and qyvalstotminus[i]==qyvalstotminus[j]:
				Uminus[2*i][2*j]=1
				Uminus[2*i+1][2*j+1]=1


			if qvalslayerminus[i]==1 and qvalslayerminus[j]==3 and qxvalstotminus[i]==qxvalstotminus[j] and qyvalstotminus[i]==qyvalstotminus[j]:
				Uminus[2*i][2*j]=1
				Uminus[2*i+1][2*j+1]=1
			if qvalslayerminus[i]==3 and qvalslayerminus[j]==1 and qxvalstotminus[i]==qxvalstotminus[j] and qyvalstotminus[i]==qyvalstotminus[j]:
				Uminus[2*i][2*j]=1
				Uminus[2*i+1][2*j+1]=1


	blockzeros=np.full((nbands,nbands),0)
	UFull=np.block([[U,blockzeros,blockzeros,blockzeros],[blockzeros,Uminus,blockzeros,blockzeros],[blockzeros,blockzeros,U,blockzeros],[blockzeros,blockzeros,blockzeros,Uminus]])



	return UFull
def MZSpin():
	nbands=np.size(qxvalstot)*2
	U=np.full((nbands,nbands),0.+0.j)
	Uminus=np.full((nbands,nbands),0.+0.j)
	for i in range(np.size(qxvalstot)):
		for j in range(np.size(qxvalstot)):

			if qvalslayer[i]==2 and qvalslayer[j]==2 and qxvalstot[i]==qxvalstot[j] and qyvalstot[i]==qyvalstot[j]:
				U[2*i][2*j]=1
				U[2*i+1][2*j+1]=1


			if qvalslayer[i]==1 and qvalslayer[j]==3 and qxvalstot[i]==qxvalstot[j] and qyvalstot[i]==qyvalstot[j]:
				U[2*i][2*j]=1
				U[2*i+1][2*j+1]=1
			if qvalslayer[i]==3 and qvalslayer[j]==1 and qxvalstot[i]==qxvalstot[j] and qyvalstot[i]==qyvalstot[j]:
				U[2*i][2*j]=1
				U[2*i+1][2*j+1]=1

	for i in range(np.size(qxvalstotminus)):
		for j in range(np.size(qxvalstotminus)):


			if qvalslayerminus[i]==2 and qvalslayerminus[j]==2 and qxvalstotminus[i]==qxvalstotminus[j] and qyvalstotminus[i]==qyvalstotminus[j]:
				Uminus[2*i][2*j]=1
				Uminus[2*i+1][2*j+1]=1


			if qvalslayerminus[i]==1 and qvalslayerminus[j]==3 and qxvalstotminus[i]==qxvalstotminus[j] and qyvalstotminus[i]==qyvalstotminus[j]:
				Uminus[2*i][2*j]=1
				Uminus[2*i+1][2*j+1]=1
			if qvalslayerminus[i]==3 and qvalslayerminus[j]==1 and qxvalstotminus[i]==qxvalstotminus[j] and qyvalstotminus[i]==qyvalstotminus[j]:
				Uminus[2*i][2*j]=1
				Uminus[2*i+1][2*j+1]=1


	blockzeros=np.full((nbands,nbands),0)
	UFull=np.block([[U,blockzeros],[blockzeros,Uminus]])



	return UFull
def EigProj():
	base=np.full((npts,nactive*4,nactive*4),0.+0.j)

	blockivc=np.full((nbands,nbands),0.+0.j)
	mztest=MZmatrix

	for i in range(npts):

		zerblockplay=np.full((nactive,nactive),0.+0.j)
		eigblock=np.block([[-1,0,0,0],[0,1,0,0],[0,0,1,0],[0,0,0,-1]])
		base[i]=np.block([[eigblock,zerblockplay,zerblockplay,zerblockplay],[zerblockplay,eigblock,zerblockplay,zerblockplay],[zerblockplay,zerblockplay,eigblock,zerblockplay],[zerblockplay,zerblockplay,zerblockplay,eigblock]])


	return base
def EigProjSpin():
	base=np.full((npts,nactive*2,nactive*2),0.+0.j)

	blockivc=np.full((nbands,nbands),0.+0.j)
	mztest=MZmatrix

	for i in range(npts):

		zerblockplay=np.full((nactive,nactive),0.+0.j)
		eigblock=np.block([[-1,0,0,0],[0,1,0,0],[0,0,1,0],[0,0,0,-1]])
		base[i]=np.block([[eigblock,zerblockplay],[zerblockplay,eigblock]])


	return base
def MZProj():
	base=np.full((npts,nactive*4,nactive*4),0.+0.j)

	blockivc=np.full((nbands,nbands),0.+0.j)
	mztest=MZmatrix

	for i in range(npts):

		U=HtotalEigVecs[i]
		base[i]=np.transpose(np.conjugate(U)).dot(mztest).dot(U)


	return base
def MZProjSpin():
	base=np.full((npts,nactive*2,nactive*2),0.+0.j)

	blockivc=np.full((nbands,nbands),0.+0.j)
	mztest=MZSpin()

	for i in range(npts):

		U=HtotalEigVecs[i]
		base[i]=np.transpose(np.conjugate(U[0:2*nbands,0:2*nactive])).dot(mztest).dot(U[0:2*nbands,0:2*nactive])


	return base
def C2ZTProj():
	base=np.full((npts,nactive*4,nactive*4),0.+0.j)

	blockivc=np.full((nbands,nbands),0.+0.j)
	mztest=MZmatrix

	for i in range(npts):

		U=HtotalEigVecs[i]
		base[i]=np.transpose(np.conjugate(U)).dot(C2zSpin.dot(T0matrix)).dot(np.conjugate(U))


	return base
def eigsystemh(H):
    w, v = np.linalg.eigh(H)
    idx = np.argsort(w)
    w = w[idx]
    v = v[:,idx]
    # print(np.abs(v[:,0]))
    return w,v
def eigsystemh(H):
    w, v = np.linalg.eig(H)
    idx = np.argsort(w)
    w = w[idx]
    v = v[:,idx]
    return w,v
# Unprojecting
def UnProjected(P):
    PUnProjected=np.full((npts,4*nbands,4*nbands),0.+0.j)

    for i in range(npts):
        U=HtotalEigVecs[i]
        PUnProjected[i]=np.conjugate(U).dot(P[i]).dot(np.transpose(U))
    return PUnProjected
def UnProjectedH(P):
    PUnProjected=np.full((npts,4*nbands,4*nbands),0.+0.j)

    for i in range(npts):
        U=HtotalEigVecs[i]
        PUnProjected[i]=U.dot(P[i]).dot(np.transpose(np.conjugate(U)))
    return PUnProjected
# Nonconstant mu
def FermiDirac(epsilon0):
	if np.abs(epsilon0)<10**(-2):
		return 1/4*(1-np.sign(epsilon0))
	return 1/2*(1-np.sign(epsilon0))
def nfunction(eps,mu):
	nelectrons=0
	for i in range(npts):
		disp=eigsystemh(eps[i])[0]

		sumtest=0

		for a in range(4*nactive):
		    nelectrons=nelectrons+1/npts*FermiDirac(disp[a]-mu)
		    sumtest=sumtest+FermiDirac(disp[a]-mu)


	return nelectrons-8
def fixChemicalPotential(eps,n0):
	mu=0
	mumin=-400
	mumax=400
	n=nfunction(eps,mu)
	i=0
	while np.abs(n-n0)>.0001:
		i=i+1
		

		if n>n0:
		    mu=mu-(mumax-mumin)/(2*2**i)
		if n<n0:
		    mu=mu+(mumax-mumin)/(2*2**i)
		n=nfunction(eps,mu)
	n=nfunction(eps,mu)
	print('n: ',n)
	return mu
##################################################################################################################################
#   Functional    #
###################

# @jit(nopython=True)
def makeFock1(FockTermId,FockTermLowerBands):
	for i in range(npts):
		for j in range(npts):
			time1=time.time()


			if np.abs(kxvals[kqindex[i,j]])<10**(-3)*ktheta and np.abs(kyvals[kqindex[i,j]])<10**(-3)*ktheta:
					continue


							

			qx1=kxvals[j]
			qy1=kyvals[j]

			if np.abs(qx1/ktheta)<10**(-4) and np.abs(qy1/ktheta)<10**(-4):
				coulombPotential=eSqOverEpsilon0*1/(2*epsilon)*dscreen*2
			else:

				modq=np.sqrt(qx1**2+qy1**2)
				coulombPotential=eSqOverEpsilon0*1/(2*epsilon)*1/modq*(1-np.exp(-dscreen*modq*2))


			
			if kqbvector[i,j]==-1:


				FockTermId[i]=FockTermId[i]- Prefactor*coulombPotential*InnerProducts[i][kqindex[i,j]].dot(np.transpose(orderId[kqindex[i,j]])).dot(np.conjugate(np.transpose(InnerProducts[i][kqindex[i,j]])))
				FockTermLowerBands[i]=FockTermLowerBands[i]-  Prefactor*coulombPotential*InnerProducts[i][kqindex[i,j]].dot(np.transpose(orderLowerBands[kqindex[i,j]])).dot(np.conjugate(np.transpose(InnerProducts[i][kqindex[i,j]])))
			else:

				if (np.abs(kxvals[kqindex[i,j]]+b1x)<10**(-3)*ktheta or np.abs(kyvals[kqindex[i,j]]-kxvals[kqindex[i,j]]*q2y/q2x-ktheta)<10**(-3) or np.abs(kyvals[kqindex[i,j]]+kxvals[kqindex[i,j]]*q2y/q2x-ktheta)<10**(-3)):


					FockTermId[i]=FockTermId[i]- Prefactor*coulombPotential*InnerProductsShell1[kqbvector[i,j]][i][kqindex[i,j]].dot(np.transpose(orderIdShell1[kqbvector[i,j]][kqindex[i,j]])).dot(np.conjugate(np.transpose(InnerProductsShell1[kqbvector[i,j]][i][kqindex[i,j]])))

					FockTermLowerBands[i]=FockTermLowerBands[i]-  Prefactor*coulombPotential*InnerProductsShell1[kqbvector[i,j]][i][kqindex[i,j]].dot(np.transpose(orderLowerBandsShell1[kqbvector[i,j]][kqindex[i,j]])).dot(np.conjugate(np.transpose(InnerProductsShell1[kqbvector[i,j]][i][kqindex[i,j]])))



				else:


					FockTermId[i]=FockTermId[i]- Prefactor*coulombPotential*InnerProductsShell1[kqbvector[i,j]][i][kqindex[i,j]].dot(np.transpose(orderId[kqindex[i,j]])).dot(np.conjugate(np.transpose(InnerProductsShell1[kqbvector[i,j]][i][kqindex[i,j]])))
					FockTermLowerBands[i]=FockTermLowerBands[i]-  Prefactor*coulombPotential*InnerProductsShell1[kqbvector[i,j]][i][kqindex[i,j]].dot(np.transpose(orderLowerBands[kqindex[i,j]])).dot(np.conjugate(np.transpose(InnerProductsShell1[kqbvector[i,j]][i][kqindex[i,j]])))
# @jit(nopython=True)
def makeFock2(FockTermId,FockTermLowerBands):
	for i in range(npts):
		for j in range(npts):
			if np.abs(kxvals[kqindex[i,j]])<10**(-3)*ktheta and np.abs(kyvals[kqindex[i,j]])<10**(-3)*ktheta:
					continue

			

			for k in range(shell1size):

				if np.sqrt((shell1momentax[k]+kxvals[j])**2+(shell1momentay[k]+kyvals[j])**2)/ktheta >2-.01 and (np.abs(kxvals[j]+b1x)<10**(-3)*ktheta or np.abs(kyvals[j]-kxvals[j]*q2y/q2x-ktheta)<10**(-3) or np.abs(kyvals[j]+kxvals[j]*q2y/q2x-ktheta)<10**(-3)):
					continue

				gx=shell1momentax[k]+shell1momentax[kqbvector[i,j]]
				gy=shell1momentay[k]+shell1momentay[kqbvector[i,j]]

				if kqbvector[i,j]==-1:
					gx=shell1momentax[k]
					gy=shell1momentay[k]

				qx1=kxvals[j]+shell1momentax[k]
				qy1=kyvals[j]+shell1momentay[k]
				modq=np.sqrt(qx1**2+qy1**2)
				coulombPotential=eSqOverEpsilon0*1/(2*epsilon)*1/modq*(1-np.exp(-dscreen*modq*2))

										
				
				xshell1=np.where(shell1momentaxround==round(gx/ktheta, 3))[0]
				yshell1=np.where(shell1momentayround==round(gy/ktheta, 3))[0]

				indexshell1=[]

				for count1 in range(len(xshell1)):
					for count2 in range(len(yshell1)):
						if xshell1[count1]==yshell1[count2]:
							indexshell1.append(xshell1[count1])
				

				

				xshell2=np.where(shell2momentaxround==round(gx/ktheta, 3))[0]
				yshell2=np.where(shell2momentayround==round(gy/ktheta, 3))[0]


				indexshell2=[]

				for count1 in range(len(xshell2)):
					for count2 in range(len(yshell2)):
						if xshell2[count1]==yshell2[count2]:
							indexshell2.append(xshell2[count1])



				if len(indexshell1)>0:

					index=indexshell1[0]

					if (np.abs(kxvals[kqindex[i,j]]+b1x)<10**(-3)*ktheta or np.abs(kyvals[kqindex[i,j]]-kxvals[kqindex[i,j]]*q2y/q2x-ktheta)<10**(-3) or np.abs(kyvals[kqindex[i,j]]+kxvals[kqindex[i,j]]*q2y/q2x-ktheta)<10**(-3)):

						FockTermId[i]=FockTermId[i]- Prefactor*coulombPotential*InnerProductsShell1[index][i][kqindex[i,j]].dot(np.transpose(orderIdShell1[index][kqindex[i,j]])).dot(np.conjugate(np.transpose(InnerProductsShell1[index][i][kqindex[i,j]])))


						FockTermLowerBands[i]=FockTermLowerBands[i]- Prefactor*coulombPotential*InnerProductsShell1[index][i][kqindex[i,j]].dot(np.transpose(orderLowerBandsShell1[index][kqindex[i,j]])).dot(np.conjugate(np.transpose(InnerProductsShell1[index][i][kqindex[i,j]])))

					else:

						FockTermId[i]=FockTermId[i]- Prefactor*coulombPotential*InnerProductsShell1[index][i][kqindex[i,j]].dot(np.transpose(orderId[kqindex[i,j]])).dot(np.conjugate(np.transpose(InnerProductsShell1[index][i][kqindex[i,j]])))


						FockTermLowerBands[i]=FockTermLowerBands[i]- Prefactor*coulombPotential*InnerProductsShell1[index][i][kqindex[i,j]].dot(np.transpose(orderLowerBands[kqindex[i,j]])).dot(np.conjugate(np.transpose(InnerProductsShell1[index][i][kqindex[i,j]])))

					
				elif len(indexshell2)>0:
					# print('hi')
					index=indexshell2[0]

					if (np.abs(kxvals[kqindex[i,j]]+b1x)<10**(-3)*ktheta or np.abs(kyvals[kqindex[i,j]]-kxvals[kqindex[i,j]]*q2y/q2x-ktheta)<10**(-3) or np.abs(kyvals[kqindex[i,j]]+kxvals[kqindex[i,j]]*q2y/q2x-ktheta)<10**(-3)):

						FockTermId[i]=FockTermId[i]- Prefactor*coulombPotential*InnerProductsShell2[index][i][kqindex[i,j]].dot(np.transpose(orderIdShell2[index][kqindex[i,j]])).dot(np.conjugate(np.transpose(InnerProductsShell2[index][i][kqindex[i,j]])))



						FockTermLowerBands[i]=FockTermLowerBands[i]- Prefactor*coulombPotential*InnerProductsShell2[index][i][kqindex[i,j]].dot(np.transpose(orderLowerBandsShell2[index][kqindex[i,j]])).dot(np.conjugate(np.transpose(InnerProductsShell2[index][i][kqindex[i,j]])))

					else:
						FockTermId[i]=FockTermId[i]- Prefactor*coulombPotential*InnerProductsShell2[index][i][kqindex[i,j]].dot(np.transpose(orderId[kqindex[i,j]])).dot(np.conjugate(np.transpose(InnerProductsShell2[index][i][kqindex[i,j]])))


						FockTermLowerBands[i]=FockTermLowerBands[i]- Prefactor*coulombPotential*InnerProductsShell2[index][i][kqindex[i,j]].dot(np.transpose(orderLowerBands[kqindex[i,j]])).dot(np.conjugate(np.transpose(InnerProductsShell2[index][i][kqindex[i,j]])))

				elif np.abs(gx/ktheta)<10**(-3) and np.abs(gy/ktheta)<10**(-3):

					FockTermId[i]=FockTermId[i]- Prefactor*coulombPotential*InnerProducts[i][kqindex[i,j]].dot(np.transpose(orderId[kqindex[i,j]])).dot(np.conjugate(np.transpose(InnerProducts[i][kqindex[i,j]])))



					FockTermLowerBands[i]=FockTermLowerBands[i]- Prefactor*coulombPotential*InnerProducts[i][kqindex[i,j]].dot(np.transpose(orderLowerBands[kqindex[i,j]])).dot(np.conjugate(np.transpose(InnerProducts[i][kqindex[i,j]])))
				else:
					print("HOW???")

# @jit(nopython=True)
def makeHartree2(HartreeTermId,HartreeTermLowerBands):

	for i in range(npts):
		for j in range(npts):

			if np.abs(kxvals[j])<10**(-3)*ktheta and np.abs(kyvals[j])<10**(-3)*ktheta:
				continue

			if np.abs(kxvals[j]+b1x)<10**(-3)*ktheta or np.abs(kyvals[j]-kxvals[j]*q2y/q2x-ktheta)<10**(-3) or np.abs(kyvals[j]+kxvals[j]*q2y/q2x-ktheta)<10**(-3):
					continue


								

			for k in range(shell1size):

				qx1=shell1momentax[k]
				qy1=shell1momentay[k]
				modq=np.sqrt(qx1**2+qy1**2)
				coulombPotential=eSqOverEpsilon0*1/(2*epsilon)*1/modq*(1-np.exp(-dscreen*modq*2))



				EdgeInnerProduct1=np.transpose(np.conjugate(HtotalEigVecs[i])).dot(HtotalEigVecsShell1[k][i])
				EdgeInnerProduct2=np.transpose(np.conjugate(HtotalEigVecs[j])).dot(HtotalEigVecsShell1[k][j])

				if i==0:
					EdgeInnerProduct1=np.transpose(np.conjugate(HtotalEigVecs[i])).dot(Kvec[k])

				if i==nptsd-1:
					EdgeInnerProduct1=np.transpose(np.conjugate(HtotalEigVecs[i])).dot(Kprimevec[k])


				if np.max(np.abs(EdgeInnerProduct1))>.98 or np.max(np.abs(EdgeInnerProduct1))>.98:
					print(i,j,EdgeInnerProduct1,EdgeInnerProduct2)


				HartreeTermLowerBands[i]=HartreeTermLowerBands[i]+1/2*(Prefactor*coulombPotential*EdgeInnerProduct1*np.trace(np.transpose(orderLowerBands[j]).dot(np.conjugate(np.transpose(EdgeInnerProduct2)))))+np.conjugate(np.transpose(1/2*(Prefactor*coulombPotential*EdgeInnerProduct1*np.trace(np.transpose(orderLowerBands[j]).dot(np.conjugate(np.transpose(EdgeInnerProduct2)))))))

				
				HartreeTermId[i]=HartreeTermId[i]+1/2*(Prefactor*coulombPotential*EdgeInnerProduct1*np.trace(np.transpose(orderId[j]).dot(np.conjugate(np.transpose(EdgeInnerProduct2)))))+np.conjugate(np.transpose(1/2*(Prefactor*coulombPotential*EdgeInnerProduct1*np.trace(np.transpose(orderId[j]).dot(np.conjugate(np.transpose(EdgeInnerProduct2)))))))
###############################################################################################################################
							

###############################################################################################################################

# Load eigenvectors and array of kxvals and kyvals - must be compatible with nptsd in save_vectors.py
# REPLACE REFERS TO MY DIRECTORIES, REPLACE WITH OWN PATHS IF USING
HtotalEigVecs=np.load('/n/holylfs/LABS/sachdev_lab/vectors2/HtotalEigVecsrepo_%d.npz'%(int(sys.argv[1]),))
HtotalEigVecsShell1=np.load('/n/holylfs/LABS/sachdev_lab/vectors2/HtotalEigVecsShell1repo_%d.npz'%(int(sys.argv[1]),))
HtotalEigVecsShell2=np.load('/n/holylfs/LABS/sachdev_lab/vectors2/HtotalEigVecsShell2repo_%d.npz'%(int(sys.argv[1]),))
HtotalEigVals=np.load('/n/holylfs/LABS/sachdev_lab/vectors2/HtotalEigValsrepo_%d.npz'%(int(sys.argv[1]),))

HtotalEigVals = HtotalEigVals.f.arr_0
HtotalEigVecs = HtotalEigVecs.f.arr_0
HtotalEigVecsShell1 = HtotalEigVecsShell1.f.arr_0
HtotalEigVecsShell2 = HtotalEigVecsShell2.f.arr_0



InnerProducts=np.full((npts,npts,4*nactive,4*nactive),0.+0.j)
InnerProductsShell1=np.full((np.size(shell1momentax),npts,npts,4*nactive,4*nactive),0.+0.j)
InnerProductsShell2=np.full((np.size(shell2momentax),npts,npts,4*nactive,4*nactive),0.+0.j)

# @jit(nopython=True)
def makeInnerProducts(InnerProducts,InnerProductsShell1,InnerProductsShell2):
	for i in range(npts):
		print('i0: ',i)
		for j in range(npts):
			InnerProducts[i][j]=np.transpose(np.conjugate(HtotalEigVecs[i])).dot(HtotalEigVecs[j])
			# print('Did1')
			for k in range(np.size(shell1momentax)):
				InnerProductsShell1[k][i][j]=np.transpose(np.conjugate(HtotalEigVecs[i])).dot(HtotalEigVecsShell1[k][j])
			# print('Did2')
			for k in range(np.size(shell2momentax)):
				InnerProductsShell2[k][i][j]=np.transpose(np.conjugate(HtotalEigVecs[i])).dot(HtotalEigVecsShell2[k][j])
			# print('Did3')

print('starting')
time0=time.time()


kxvals=np.load('kxvalsrepo.npy')
kyvals=np.load('kyvalsrepo.npy')
kqindex=np.load('kqindexrepo.npy')
kqbvector=np.load('kqbvectorrepo.npy')
kminusindex=np.load('kminusindexrepo.npy')
kc3index=np.load('kc3indexrepo.npy')
kc2yindex=np.load('kc2yindexrepo.npy')


makeInnerProducts(InnerProducts,InnerProductsShell1,InnerProductsShell2)

C3matrix=C3()
C2zmatrix=C2z()
C3Spinmatrix=C3Spin()
C2zSpinmatrix=C2zSpin()
T0matrix=T0()
T0Spinmatrix=T0Spin()
U1vmatrix=U1v()
Chiralmatrix=Chiral()
Spinmatrix=Spin()
MZmatrix=MZ()

C3C2product=C3matrix.dot(C2zmatrix)
C3C3C2product=C3matrix.dot(C3matrix).dot(C2zmatrix)
C3C3product=C3matrix.dot(C3matrix)
C3T0product=C3matrix.dot(T0matrix)
C3C3T0product=C3matrix.dot(C3matrix).dot(T0matrix)

print('NAN?: ',np.isnan(np.sum(HtotalEigVecsShell1)))

# Arrays to store quantities useful later

prevELowerBands=np.full((npts,4*nactive,4*nactive),0.+0.j)

KineticTerm=[]
# Make kinetic term from eigenvalues of continuum model
for i in range(npts):
    KineticTerm.append(np.diag(HtotalEigVals[i]))



print('Initial Orders Go')
orderLowerBands=LowerBands()



print('Start Orders Go')


orderId=Identity()
orderIdShell1=Identityshell1()
orderIdShell2=Identityshell2()
print('Shell 1 Orders Go')


time0=time.time()
print('time: ',time.time()-time0)
print('Shell2 Orders Go')

print('Shell2 Part 1')
time0=time.time()
print('time: ',time.time()-time0)
print('Shell2 Part 2')
orderLowerBandsShell1=MakeShell1PC2z(orderLowerBands)
orderLowerBandsShell2=MakeShell2PC2z(orderLowerBands)

Kvec=np.full((np.size(shell1momentax),4*nbands,4*nactive),0.+0.j)
Kprimevec=np.full((np.size(shell1momentax),4*nbands,4*nactive),0.+0.j)

for i in range(np.size(shell1momentax)):
	print('shell: ',i)

	HtotalEigVecsShell1x=[]
	for j in range(npts):
		# print('pt: ',j)


		unew=np.full(nbands,0.+0.j)
		hplusnew=np.full((nbands,nactive),0.+0.j)
		hminusnew=np.full((nbands,nactive),0.+0.j)



		for k in range(np.size(qvalslayer)):

			checkplusx=np.where(np.round(qxvalstot/ktheta,3)==np.round((qxvalstot[k]-shell1momentax[i])/ktheta,3))
			checkplusy=np.where(np.round(qyvalstot/ktheta,3)==np.round((qyvalstot[k]-shell1momentay[i])/ktheta,3))

			checkpluslayer=np.where(np.round(qvalslayer,3)==np.round(qvalslayer[k],3))
			intersectplus=np.intersect1d(checkpluslayer,np.intersect1d(checkplusx, checkplusy))


			if np.size(intersectplus)>0:

				if j==0:

					Kvec[i][2*k,:]=HtotalEigVecs[j][2*intersectplus[0],:]
					Kvec[i][2*k+1,:]=HtotalEigVecs[j][2*intersectplus[0]+1,:]
					Kvec[i][2*nbands+2*k,:]=HtotalEigVecs[j][2*nbands+2*intersectplus[0],:]
					Kvec[i][2*nbands+2*k+1,:]=HtotalEigVecs[j][2*nbands+2*intersectplus[0]+1,:]

				if j==nptsd-1:

					Kprimevec[i][2*k,:]=HtotalEigVecs[j][2*intersectplus[0],:]
					Kprimevec[i][2*k+1,:]=HtotalEigVecs[j][2*intersectplus[0]+1,:]
					Kprimevec[i][2*nbands+2*k,:]=HtotalEigVecs[j][2*nbands+2*intersectplus[0],:]
					Kprimevec[i][2*nbands+2*k+1,:]=HtotalEigVecs[j][2*nbands+2*intersectplus[0]+1,:]

				

			else:


				if j==0:

					Kvec[i][2*k,:]=0
					Kvec[i][2*k+1,:]=0
				if j==nptsd-1:

					Kprimevec[i][2*k,:]=0
					Kprimevec[i][2*k+1,:]=0


			checkminusx=np.where(np.round(qxvalstotminus/ktheta,3)==np.round((qxvalstotminus[k]-shell1momentax[i])/ktheta,3))
			checkminusy=np.where(np.round(qyvalstotminus/ktheta,3)==np.round((qyvalstotminus[k]-shell1momentay[i])/ktheta,3))

			checkminuslayer=np.where(np.round(qvalslayerminus,3)==np.round(qvalslayerminus[k],3))
			intersectminus=np.intersect1d(checkminuslayer,np.intersect1d(checkminusx, checkminusy))

			if np.size(intersectminus)>0:


				if j==0:

					Kvec[i][nbands+2*k,:]=HtotalEigVecs[j][nbands+2*intersectminus[0],:]
					Kvec[i][nbands+2*k+1,:]=HtotalEigVecs[j][nbands+2*intersectminus[0]+1,:]
					Kvec[i][3*nbands+2*k,:]=HtotalEigVecs[j][3*nbands+2*intersectminus[0],:]
					Kvec[i][3*nbands+2*k+1,:]=HtotalEigVecs[j][3*nbands+2*intersectminus[0]+1,:]

				if j==nptsd-1:

					Kprimevec[i][nbands+2*k,:]=HtotalEigVecs[j][nbands+2*intersectminus[0],:]
					Kprimevec[i][nbands+2*k+1,:]=HtotalEigVecs[j][nbands+2*intersectminus[0]+1,:]
					Kprimevec[i][3*nbands+2*k,:]=HtotalEigVecs[j][3*nbands+2*intersectminus[0],:]
					Kprimevec[i][3*nbands+2*k+1,:]=HtotalEigVecs[j][3*nbands+2*intersectminus[0]+1,:]


			else:
				# if j==0 :
				# 	print('0')

				if j==0:

					Kvec[i][nbands+2*k,:]=0
					Kvec[i][nbands+2*k+1,:]=0
					Kvec[i][3*nbands+2*k,:]=0
					Kvec[i][3*nbands+2*k+1,:]=0

				if j==nptsd-1:

					Kprimevec[i][nbands+2*k,:]=0
					Kprimevec[i][nbands+2*k+1,:]=0
					Kprimevec[i][3*nbands+2*k,:]=0
					Kprimevec[i][3*nbands+2*k+1,:]=0

print('Finished')

for it in range(iterations):

	print('shell1')
	

	orderLowerBandsShell1=MakeShell1P(orderLowerBands)
	orderLowerBandsShell2=MakeShell2P(orderLowerBands)

	



	print('Final Symmetries')
	# print('c3: LowerBands ',C3Symmetry(UnProjected(orderLowerBands)))
	# print('c2z: LowerBands ',C2zSymmetry(UnProjected(orderLowerBands)))
	# print('t: LowerBands ',TSymmetry(UnProjected(orderLowerBands)))



	FockTermId=np.full((npts, 4*nactive,4*nactive),0.+0.j)
	HartreeTermId=np.full((npts, 4*nactive,4*nactive),0.+0.j)

	FockTermLowerBands=np.full((npts, 4*nactive,4*nactive),0.+0.j)
	HartreeTermLowerBands=np.full((npts, 4*nactive,4*nactive),0.+0.j)

	print('iteration: ',it)
	if it>0:
	    print('time of last iteration: ',time.time()-timestart)
	timestart=time.time()

	print('Starting Construction of Hartree-Fock functional')
	timestart=time.time()
	makeFock2(FockTermId,FockTermLowerBands)
	makeFock1(FockTermId,FockTermLowerBands)
	makeHartree2(HartreeTermId,HartreeTermLowerBands)
	print('time: ',time.time()-timestart)
	print('did Fock Term')
	print(np.sum(np.abs(FockTermLowerBands)))

	colors=[]
	maxmin=[]



	energyFilledbands=1/2*AMBZ*A/(N*npts*(2*np.pi)**2)*(2*np.einsum('iab,iab',KineticTerm,orderLowerBands)+np.einsum('iab,iab',FockTermLowerBands,orderLowerBands)+2*np.einsum('iab,iab',FockTermId,orderLowerBands)+np.einsum('iab,iab',HartreeTermLowerBands,orderLowerBands)+2*np.einsum('iab,iab',HartreeTermId,orderLowerBands))

	energySubtraction=1/2*AMBZ*A/(N*npts*(2*np.pi)**2)*(2*np.einsum('iab,iab',KineticTerm,orderId)+np.einsum('iab,iab',FockTermId,orderId)+2*np.einsum('iab,iab',FockTermId,orderId)+np.einsum('iab,iab',HartreeTermId,orderId)+2*np.einsum('iab,iab',HartreeTermId,orderId))


	print('LowerBands: ',energyFilledbands)

	muIdentity=np.full((npts,4*nactive,4*nactive),0.+0.j)
	for i in range(npts):
		muIdentity[i]=np.identity(4*nactive)






	print('Finished Construction of Hartree-Fock functional')



	HtotalLowerBands=KineticTerm+FockTermLowerBands+HartreeTermLowerBands

	muLowerBands=fixChemicalPotential(HtotalLowerBands,0)

	print('LowerBands: ',muLowerBands)


	HtotalLowerBands=HtotalLowerBands-muLowerBands*muIdentity



	# print('hamiltonian symmetry LowerBands: ',C3SymmetryH(UnProjectedH(HtotalLowerBands)))
	print('hamiltonian symmetry LowerBands: ',TSymmetryH(UnProjectedH(HtotalLowerBands)))
	# print('hamiltonian symmetry LowerBands: ',C2zSymmetryH(UnProjectedH(HtotalLowerBands)))





	edifLowerBands=0

	edifre=0
	edifim=0

	nLowerBands=0



	for i in range(npts):


		ULowerBandsv1s1=eigsystemh(HtotalLowerBands[i][0:nactivehalf*2,0:nactivehalf*2])[1]
		ULowerBandsv2s1=eigsystemh(HtotalLowerBands[i][nactivehalf*2:2*nactivehalf*2,nactivehalf*2:2*nactivehalf*2])[1]
		ULowerBandsv1s2=eigsystemh(HtotalLowerBands[i][2*nactivehalf*2:3*nactivehalf*2,2*nactivehalf*2:3*nactivehalf*2])[1]
		ULowerBandsv2s2=eigsystemh(HtotalLowerBands[i][3*nactivehalf*2:4*nactivehalf*2,3*nactivehalf*2:4*nactivehalf*2])[1]

		Pblock=np.full((nactivehalf*2,nactivehalf*2),0.+0.j)

		v=np.block([[ULowerBandsv1s1,Pblock,Pblock,Pblock],[Pblock,ULowerBandsv2s1,Pblock,Pblock],[Pblock,Pblock,ULowerBandsv1s2,Pblock],[Pblock,Pblock,Pblock,ULowerBandsv2s2]])
		w=np.hstack((eigsystemh(HtotalLowerBands[i][0:nactivehalf*2,0:nactivehalf*2])[0],eigsystemh(HtotalLowerBands[i][nactivehalf*2:2*nactivehalf*2,nactivehalf*2:2*nactivehalf*2])[0],eigsystemh(HtotalLowerBands[i][2*nactivehalf*2:3*nactivehalf*2,2*nactivehalf*2:3*nactivehalf*2])[0],eigsystemh(HtotalLowerBands[i][3*nactivehalf*2:4*nactivehalf*2,3*nactivehalf*2:4*nactivehalf*2])[0]))
		idx = np.argsort(w)
		w = w[idx]
		ULowerBands = v[:,idx]
		wprime=1/2*(np.abs(np.full((np.size(w)),1.)-np.sign(w)))
		Chiprime=np.diag(wprime)


		nLowerBands=nLowerBands+1/npts*np.trace(Chiprime)

		orderLowerBands[i]=np.conjugate(ULowerBands).dot(Chiprime).dot(np.transpose(ULowerBands))


		for j in range(np.size(w)):
			if np.abs(w[j])<5*10**(-1):
				wprime[j]=1/2
		Chiprime=np.diag(wprime)


		orderLowerBands[i]=np.conjugate(ULowerBands).dot(Chiprime).dot(np.transpose(ULowerBands))









		edifLowerBands=edifLowerBands+np.abs(np.sum(prevELowerBands[i]-orderLowerBands[i]))
		prevELowerBands[i]=orderLowerBands[i]
	np.save('/n/home11/mchristos/TTG/epsvary/eps10noPHtest2HtotalLowerBandsrepo_%d'%(int(sys.argv[1])),HtotalLowerBands)

	print('energy difference LowerBands: ',edifLowerBands)


